text:
module TIFFtoRegression

############################################################################
# Reading the values of the TIFF according to TIFF v0.6

function FieldValues(ValueOffset::Int64,FileNumber::Array{Int16},FieldCount::Int64,FieldType=4::Int64) #v2
# =========================== Value or Pointer ? ===============================
ValuePointer=0::Int64
if FieldCount==1 # Simplified criterion of "can hold that much info on 4 bytes"
  Value=ValueOffset
else
  ValuePointer=ValueOffset+1 # correction due to the first being 0 and not 1
end

if ValuePointer==0
  DesiredV=Value
else
  #println("$ValuePointer, Count=$FieldCount, Type=$FieldType")
  DesiredV=zeros(1,FieldCount)
  for i=1:FieldCount
    Vi=0
    if FieldType == 2 || FieldType==1 #ASCII or BYTE
      FieldTypeByte=1
    elseif FieldType==3 # Short
      FieldTypeByte=2
    elseif FieldType==4 # Long
      FieldTypeByte=4
    elseif FieldType==5 # Rational
      FieldTypeByte=8
    end

    for jj=1:FieldTypeByte
      Vjj=FileNumber[ValuePointer+FieldTypeByte*(i-1)+jj-1]*256^(jj-1)
      Vi=Vi+Vjj
    end
    DesiredV[i]=convert(Int64,Vi)
  end
end

return DesiredV
end
#############################################################################

# Reading a TIFF (also able to read special Dantec TIFF that are pseudo Std)

function AquiTiff(FilePath::String)
  # Opening File
  FileString=read(FilePath)
  L=length(FileString)
  FileNumber=zeros(Int16,L,1)
  for i=1:L
    NumberI=convert(Int16,FileString[i])
    FileNumber[i]=NumberI
  end
  # ============================ Reading Header ================================
  IFD=FileNumber[5]+1+FileNumber[6]*256+FileNumber[7]*256^2+FileNumber[8]*256^3 # indicates where the Image File Directory starts (i e where the meta data start)
  # +1 in the indexes is a consequence of the starting count from 0 to 1 in tiff ==> JuliaVector
  if FileNumber[1]==73 & FileNumber[2]==73 # Small endian (determines the order of the bites from bigger to smaller)
  # ======================== Reading Metadata: ================================
  # ------------------------- Opening Values ----------------------------------
    Entries=FileNumber[IFD]+256*FileNumber[IFD+1] # Number of Entries in the IFD
    # Opening Values that might be in the Metadata
    Height=0::Int64
    Width=0::Int64
    BPS=0::Int64 #Bit per Sample
    CompType=0::Int64 #Type of Compression
    BlackWhite=0::Int64 #Photometric Interpretation # 0=White is Zero; 1<=Black is Zero
    StripOffsets=0::Int64 #Places where data is stored
    NFragData=0::Int64 # Number of those Places
    Orientation=0::Int64 # 1=Start Top-left
    SPPx=0::Int64 # Sample per Pixel
    RowsPStrip=0::Int64 # Rows per data strip: ceilling(Height/RowsPStrip)=Number of strip (p19)
    StripByteCount=0::Int64 #For each strip the number of bytes in that strip AFTER compression
    XRes=0::Int64 # X-Resolution : nbr of px per Resolution unit (cm/inch see 296)
    YRes=0::Int64 # Y-Resolution
    PlanConfig=0::Int64 # Planar Configuration ???
    ResUnit=0::Int64 # Resolution Unit : abstract (=0), cm (=1) or inch (=2)
    Software=0::Int64 # Software used to create the Tiff
    SBCType=0::Int64
    StripOffsets=[]
    NBytesPPx=0::Int64 # Spec Dantec Dynamics
    # --------------------------------------------------------------------------

    for i=1:Entries # All MetaData there are
      # ----------------------- Initiating Entry--------------------------------
      # Delimiting entry
      StartE=IFD+12*(i-1)+2
      EndE=StartE+11
      EntryN=FileNumber[StartE:EndE] #Taking one entry of 12 hexa bytes
      # Reading Entry
      FieldID=EntryN[1]+EntryN[2]*256
      FieldType=EntryN[3]+EntryN[4]*256
      FieldCount=EntryN[5]+EntryN[6]*256+EntryN[7]*65536+EntryN[8]*16777216
      ValueOffset=EntryN[9]+EntryN[10]*256+EntryN[11]*65536+EntryN[12]*16777216
      # Is either the value of the entry or a pointer towards the place where the values are kept if the entry can't be held on 8bits
      # ----------------------- Switch on FieldID ------------------------------
      if FieldID==254 # NewSubFileType
        NewSubFileType=FieldValues(ValueOffset,FileNumber,FieldCount)
        # 2 : Multipage
        # 4 : Defines a transparency mask
      elseif FieldID==256 # Image Width= n° of column in the Image
        Width=FieldValues(ValueOffset,FileNumber,FieldCount)
      elseif FieldID==257 # Height
        Height=FieldValues(ValueOffset,FileNumber,FieldCount)
      elseif FieldID==258 # Bits/sample ==> number of shade of gray allowable
        BPS=FieldValues(ValueOffset,FileNumber,FieldCount)
      elseif FieldID==259 # Compression Type
        CompType=FieldValues(ValueOffset,FileNumber,FieldCount)
      elseif FieldID==262 # Photometric Interpretation
        # 0=White is Zero; 1<=Black is Zero
        BlackXhite=FieldValues(ValueOffset,FileNumber,FieldCount)
      elseif FieldID==273 # StripOffsets= Position of all data StripOffsets
        NFragData=FieldCount
        StripOffsets=FieldValues(ValueOffset,FileNumber,FieldCount)
      elseif FieldID==274 # Orientation of the Picture
        # 1=Start Top-left
        Orientation=FieldValues(ValueOffset,FileNumber,FieldCount)
      elseif FieldID==277 # Sample/px (3 or more ==> Number of colors allowable)
        SPPx=FieldValues(ValueOffset,FileNumber,FieldCount)
      elseif FieldID==278 # Rows per Strip
        RowsPStrip=FieldValues(ValueOffset,FileNumber,FieldCount)
      elseif FieldID==279 # StripByteCount
        StripByteCount=FieldValues(ValueOffset,FileNumber,FieldCount)
        SBCType=FieldType
      elseif FieldID==282 # X-Resolution
        XRes=FieldValues(ValueOffset,FileNumber,FieldCount)
      elseif FieldID==283 # Y-Resolution
        YRes=FieldValues(ValueOffset,FileNumber,FieldCount)
      elseif FieldID==284 # Planar Configuration
        PlanConfig=FieldValues(ValueOffset,FileNumber,FieldCount)
      elseif FieldID==285
        PageName=FieldValues(ValueOffset,FileNumber,FieldCount,FieldType)
      elseif FieldID==296 # Resolution Unit (cm/inch)
        ResUnit=FieldValues(ValueOffset,FileNumber,FieldCount)
      elseif FieldID==297
        PageNumber=FieldValues(ValueOffset,FileNumber,FieldCount,FieldType)
      elseif FieldID==305 # Software
        #Software=FieldValues(ValueOffset,FileNumber,FieldCount)
      elseif FieldID==32997
        #Specific to DynamicStudio, Dantec
        # Number of bytes/ px? TBC
        NBytesPPx=FieldValues(ValueOffset,FileNumber,FieldCount)
        if SPPx==0
          SPPx=2 #NBytesPPx/8 # Works?
          SPPx=convert(Int64,SPPx)
        end
      end  # Metadata Aquisition
    end
    EndIFD=IFD+12*(Entries-1)+13
  elseif FileNumber[1]==77 & FileNumber[2]==77 #big endian
    println("Big Endian Still to be implemented, nearly same as small endian")
  end
  # ======================= Creating the Img Matrix ============================
  # ----------------------- Finding and Ordering Px ---------------------------
  PixelString=[]
  PixelNumber=ones(1,length(FileNumber))
  CountLeng=0
  for i=1:NFragData
    PxStrip=FileNumber[convert(Int64,StripOffsets[i]):convert(Int64,StripOffsets[i]+StripByteCount[i])]
    l=length(PxStrip)
    PixelNumber[1+1+CountLeng:1+l+CountLeng]=PxStrip[1:l]
    CountLeng=CountLeng+l
  end
  PixelNumber=PixelNumber[1:CountLeng] # cutting the exedentary ones (due to the last strip being sometimes shorter than the others)
  # ------------------------ Making Order in the Px ----------------------------
  ImgArray=zeros(Int64,Height,Width,SPPx)
  for i=1:Height
    for  j=1:Width
      for k=1:SPPx
        if NBytesPPx==0
          Pxijk=((i-1)*Width+(j-1))*SPPx+k
        else # Why ? 171 dtmed empirically for Dynamic Studio Tif)
          #Pxijk=((i-1)*Width+(j-1))*SPPx+k
          Pxijk=((i-1)*Width+(j+convert(Int64,floor(171*i/Height))-1))*SPPx+k # ???????????????????????????????????????????????????
        end
        ImgArray[i,j,k]=PixelNumber[Pxijk]
      end
    end
  end
  #--------------------- Making it Gray Matrix ---------------------------------
  if NBytesPPx==0 # Suposedly "Normal Img or Color Img"
    #GrayArray=sum(ImgArray,3)/SPPx
    GrayArray=ImgArray[:,:,3]
  else
    GrayArray=ImgArray[:,:,1]*256+ImgArray[:,:,2] # Dantec Spec
  end
return GrayArray[:,:,1]
end

###########################################################################

# Loading a Serie of TIFF from Dantec Dynamics

function AquiSerieD(Folder::String, File0::String, NbrImg=189::Int64)
  # ========================== Making Img Path =================================
  # ----------------Cutting File0 (Dantec DynamicsStudio Naming)----------------
  j=1
  while File0[j]!='.'
    j=j+1
  end
  p1=j
  #println(p1)
  BaseName=File0[1:p1-1] # Input by User exporting Img
  j=j+1
  while File0[j]!='.'
    j=j+1
  end
  p2=j
  FileSerie=File0[p1+1:p2-1] #Auto Serie Numbering
  NumberInSerie=File0[p2+1:end] # Auto Img Numbering
  N=NbrImg-1 # Numbering Starts at 0 not 1
  ImgArray=[]
  OneAdress=[]
  # ------------------------ Making Image Number String -------------------------
  for i=0:N
    # Finding k
    LNIS=length(NumberInSerie)
    k=LNIS
    while floor(i/(10.0^k))<1 && k>=0
      k=k-1
    end
    # completing the NumberIn Serie String
    NIS=""
    ik=i
    for kk=LNIS:-1:1
      if kk>k+1
        NIS="$(NIS)0"
      else
        V=convert(Int64,floor(ik/(10^(kk-1))))
        NIS="$(NIS)$V"
        ik=ik-V*10^(kk-1)
      end
    end
    OneAdress="$(Folder)\\$(BaseName).$(FileSerie).$(NIS).tif"
    #println(OneAdress)
    # ========================= Aquiring Images ==================================
    OneImg=AquiTiff(OneAdress)
    if isempty(ImgArray)
      ImgArray=zeros(Int64,size(OneImg,1),size(OneImg,2),N+1)
      ImgArray[:,:,1]=OneImg
    else
      ImgArray[:,:,i+1]=OneImg
    end
  end
return ImgArray # to assure Julia wont print the Array
end
#############################################################################

# Loading a Serie of TIFF automaticaly renamed with Windows 10

function AquiSerieW(Folder::String, File0::String, NbrImg=189::Int64)
  # ========================== Making Img Path =================================
  # ----------------Cutting File0 (Dantec DynamicsStudio Naming)----------------
  j=1
  while File0[j]!='('
    j=j+1
  end
  p1=j
  #println(p1)
  BaseName=File0[1:p1-1] # Input by User + 1 space (std windows renumerotation)
  j=j+1

  ImgArray=[]
  OneAdress=[]
  # ------------------------ Making Image Number String -------------------------
  for i=1:NbrImg
    OneAdress="$(Folder)\\$(BaseName)($i).tif"
    #println(OneAdress)
    # ========================= Aquiring Images ==================================
    OneImg=AquiTiff(OneAdress)
    if isempty(ImgArray)
      ImgArray=zeros(Int64,size(OneImg,1),size(OneImg,2),NbrImg)
      ImgArray[:,:,1]=OneImg
    else
      ImgArray[:,:,i]=OneImg
    end
  end
return ImgArray # to assure Julia wont print the Array
end

#############################################################################

# Creating a Special type to stock the start of a line and its values
type Line # Creates type Line
  YValues::Vector
  XStart::Int64
end
##############################################################################

# Detecting the Lines in an image matrix (also retourning a Binary matrix)

function DetectLine(Img::Matrix)
  # only for bicolor ImgLine
  BiMat=zeros(Bool,size(Img))
  BiMat0=zeros(Bool,size(Img))
  Max=maximum(Img)
  min=minimum(Img)
  moy=(Max+min)/2
# --------------------------BiColoring----------------------------------
  k=0
  for i=1:size(Img,1) #BiLeveling
    for j=1:size(Img,2)
      if Img[i,j]<=moy
        BiMat[i,j]=0
        BiMat0[i,j]=0
      else
        BiMat[i,j]=1
        BiMat0[i,j]=1
      end
    end
  end

SumMat=sum(BiMat)
H=size(BiMat,1)
L=size(BiMat,2)
NPx=H*L
if SumMat<NPx/2
  # 1=Line, 0=Background
  LinePx=true
  BckGrnd=false #BckGrnd LinePx
else
  # 0=Line, 1=Background
  LinePx=false
  BckGrnd=true
end
CountLinePx=abs(SumMat-NPx) # nbr of Px that are part of a line
Lines=[] # need Array to stock the XStart
Xcount=1::Int64
Split=false
Ymean2=H
YmeanM1=H
Xcount2=1
while Xcount<L+1
  x=Xcount
 # 1 Line code  !
  if Split
    y=Ymean2
    x=Xcount2
    Xcount=Xcount-1
    Split=false
  else
    y=H
    x=Xcount
  end
  ConLine=true
  # ========================= Following the Line =========================
  while y>=2 && ConLine# going up to top
    # find [Background; Line]
    YmeanV=[]
    if [LinePx; BckGrnd]==[BiMat[y-1,x]; BiMat[y,x]] # Line Start Detect
      #Xcount=Xcount-1 #Stay in the same column and check further
      # Special Start in case of Split
      if YmeanM1 != H
        push!(YmeanV,YmeanM1)
        XStart=x-1
        YmeanM1=H
      else #End Spec Split
        XStart=x
      end
      y=y-1
      YStart=y
      # Start of a Line
      while BiMat[y,x]==LinePx
        y=y-1
      end
      YEnd=y+1
      # following that one line
      Width=abs(YStart-YEnd)+1
      BiMat[YEnd:YStart,x]=repmat([BckGrnd],Width,1)
      Ymean=(YStart+YEnd)/2
      push!(YmeanV,Ymean)
      while ConLine && x<L
        x=x+1
        if YStart<H && YEnd>1
          if sum(BiMat[YEnd-1:YStart+1,x].==repmat([BckGrnd],Width+2,1))==Width+2
            # Line Stop
            ConLine=false
            YStart=0
            YEnd=0
            x=Xcount
            Xcount=Xcount-1
          else
            # Line continue
            if BiMat[YEnd-1:YStart+1,x]==[BckGrnd;repmat([LinePx],Width,1);BckGrnd]
             # Line constant in Value and Width
            else
              # Line continue but changed
              if (BiMat[YStart:YStart+1,x]==[LinePx;LinePx])||(BiMat[YEnd-1:YEnd,x]==[LinePx;LinePx])
                # Line increase or decrease
                if (BiMat[YStart:YStart+1,x]==[LinePx;LinePx]) && !(BiMat[YEnd-1:YEnd,x]==[LinePx;LinePx])
                  # Decreasing
                  YEnd=YEnd+minimum(findin(BiMat[YEnd:YStart,x].==repmat([LinePx],Width,1),true))-1
                  yy=YStart
                  while BiMat[yy+1,x]==LinePx && yy<H
                    yy=yy+1
                  end
                  YStart=yy
                elseif !(BiMat[YStart:YStart+1,x]==[LinePx;LinePx]) && (BiMat[YEnd-1:YEnd,x]==[LinePx;LinePx])# Increasing
                  YStart=YEnd+maximum(findin(BiMat[YEnd:YStart,x].==repmat([LinePx],Width,1),true))-1
                  # YEnd
                  yy=YEnd
                  while BiMat[yy-1,x]==LinePx && yy>1
                    yy=yy-1
                  end
                  YEnd=yy
                elseif (BiMat[YStart:YStart+1,x]==[LinePx;LinePx]) && (BiMat[YEnd-1:YEnd,x]==[LinePx;LinePx]) # Increse in up and down
                  if sum(BiMat[YEnd-1:YStart+1,x].==LinePx)==(YStart+1-(YEnd-1)+1)# No hole
                    # Increse in Width in both side
                    ydown=YStart
                    while BiMat[ydown+1,x]==LinePx && ydown<H
                      ydown=ydown+1
                    end
                    YStart=ydown
                    yup=YEnd
                    while BiMat[yup-1,x]==LinePx && yup>1
                      yup=yup-1
                    end
                    YEnd=yup
                  else # there is a hole
                    # Split in two curves
                    Split= true
                    PosHole=findin(BiMat[YEnd-1:YStart+1,x].==BckGrnd,true)
                    # Down Part Continues
                    YEnd=(YEnd-1+minimum(PosHole))-1 # (below hole)
                    ydown=YStart
                    while BiMat[ydown+1,x]==LinePx && ydown<H
                      ydown=ydown+1
                    end
                    YStart=ydown
                    # Up Going part
                    YStart2=(YEnd-1+maximum(PosHole))+1
                    Ymean2=YStart2+1
                    Xcount2=x
                    YmeanM1=Ymean # kept for next line
                  end
                end
              else
                #Line Width decrease within previous limits
                PxPos=findin(BiMat[YEnd:YStart,x].==repmat([LinePx],Width,1),true)
                minPxPos=minimum(PxPos)
                MaxPxPos=maximum(PxPos)
                YStart=YEnd+MaxPxPos-1
                YEnd=YEnd+minPxPos-1
              end
            end
          end
        else
          # border effect
        end
        if !ConLine
          Width=0
        else
          Width=abs(YStart-YEnd)+1
          BiMat[YEnd:YStart,x]=repmat([BckGrnd],Width,1)
          Ymean=(YEnd+YStart)/2
          push!(YmeanV,Ymean)
        end
      end #start of a line detected
    else
      y=y-1 # No Line, go on climbing
    end
    if length(YmeanV)>1 #!isempty(YmeanV)
      Line1=Line(H+1-YmeanV,XStart) # Assuming left bottom axis =0
      #then complete Lines
      push!(Lines,Line1)
      ConLine=false # end current line
    end
  end
    Xcount=Xcount+1
end
# from bottom to top
# 1 First column -> detect LinePx
# 1b if no LinePx go to the right
# 2 continue climbing until Px!=LinePx
# 3 from there moove to the right following the line
# 3a check same Px and Px+1 and Px-1
# 3b if Px+-1 check also Px+-2, etc
#    if no Px: line stop
# 3c repeat 3 until line stop
# 4 remoove selected Px
# 5 repeat operation until top right
return BiMat0, Lines
end

#############################################################################

# Allowing to plot the new type created, eventually on the matrix (thanks to PyPlot)

function PlotLine(Lines::Array{Any},H=0::Int64)
  for i=1:length(Lines)
    XS=Lines[i].XStart
    YV=Lines[i].YValues
    if H==0
      plot(XS:XS+length(YV)-1,YV,"-")
    else
      plot(XS-1:XS-1+length(YV)-1,H-YV,"-")
    end
  end
end

############################################################################

# Creating the next function for the Multiple Regression Function below

function NextFun(Xin::Vector,RegFuns="Poly"::String,i=1::Int64)
  n=length(Xin)
  if RegFuns=="Poly"
    Xout=Xin.^(i)
    XT=("Poly Order $i")
  elseif RegFuns=="Trigo"
    if isinteger(i/2)
      Xout=cos.(Xin/maximum(Xin)*(i/2))
      XT=("Standirdize cos($i x/2)")
    else
      Xout=cos.(Xin/maximum(Xin)/((i+1)/2))
      XT=("Standirdize cos($(1/(i+1)) x/2)")
    end
  elseif RegFuns=="Frac"
    F=findin(Xin,0)
    for j=1:length(F)
      X[F[j]]=mean(abs(Xin))*10^-5
    end
    Xout=1./Xin.^i
    if abs(minimum(Xout))< 2e-15 # all becoming too small
      println("Elements becoming too small X=ones")
      Xout=ones(size(Xin))
      XT=("Ones")
    else
      XT=("Frac Order $i")
    end
  elseif RegFuns=="Expo"
    Xout=exp.(i*Xin/maximum(Xin))
    XT=("Standardize exp($i x)")
  elseif RegFuns=="Loga"
    Xout=log.(abs(Xin+(i-2)*sqrt(maximum(Xin))))
    XT=("log(x+$(i-2) sqrt(Max))")
  elseif RegFuns=="Mixte"
    NFam=5 #nbr of family used (code written could be improoved)
    if isinteger((i+NFam-1)/NFam)
      k=(i+NFam-1)/NFam
      Xout=Xin.^k
      XT=("Poly Order $k") #Poly
    elseif isinteger((i+NFam-2)/NFam)
      k=(i+NFam-2)/NFam
      F=findin(Xin,0)
      for j=1:length(F)
        X[F[j]]=mean(abs(Xin))*10^-5
      end
      Xout=1./Xin.^k
      if abs(minimum(Xout))< 2e-15 # all becoming too small
        println("Elements becoming too small X=rand")
        Xout=rand(size(Xin))
      else
        XT=("Frac Order $k")
      end #Frac
    elseif isinteger((i+NFam-3)/NFam)
      k=(i+NFam-3)/NFam
      #NTrigo=4
      if isinteger(k/2)
        Xout=cos.(Xin/maximum(Xin)*(k/2))
        XT=("Standirdize cos($k x/2)")
      else
        Xout=cos.(Xin/maximum(Xin)/(k+1)/2)
        XT=("Standirdize cos($(1/(k+1)) x/2)")
      end #Trigo
    elseif isinteger((i+NFam-4)/NFam)
      k=(i+NFam-3)/NFam
      Xout=exp.(k*Xin/maximum(Xin))
      XT=("Standardize exp($k x)") #Expo
    elseif isinteger((i+NFam-5)/NFam)
      k=(i+NFam-5)/NFam
      Xout=log.(abs(Xin+(k-2)*sqrt(maximum(Xin))))
      XT=("log(x+$(k-2) sqrt(Max))")
    end
  elseif RegFuns=="Custom" # Would their be another interesting family
    println("What next fun do you want ? (Existing type: Poly, Trigo, Frac, Expo or Loga) \n Your choice is :")
    SwCu=input() #Switch Custom
    if SwCu=="Poly"
      println("Order of the wanted term : y=x^{Order} ? \n Order=")
      u=parse(Float64,input())
      if isinteger(u)
        Xu=Xin
      else
        Xu=Xin
        Xu[find(Xu.<0)]=zeros(size(find(Xu.<0)))
        println("All smaller than 0 elements have been replaced by zero for Poly Order $u")
      end
      Xout=Xu.^(u)
      XT=("Poly Order $u")
    elseif SwCu=="Trigo"
      println("Normal or Inverse ? \n Your choice is :")
      TTrig=input()
      if TTrig=="Normal"
        println("What type do you want ? (Possible choices: cos, sin, tan or cot) \n Your choice is :")
        TrigN=input()
        if TrigN=="sin"
          pritnln("What Angulat fequency (omega) do you want ? \n omega =")
          omega=parse(Float64,input())
          println("What Dephasage (phi) do you want ? \n phi =")
          phi=parse(Float64,input())
          Xout=sin.(omega*Xin+phi)
          XT=("sin($omega x+$phi)")
        elseif TrigN=="cos"
          pritnln("What Angulat fequency (omega) do you want ? \n omega =")
          omega=parse(Float64,input())
          println("What Dephasage (phi) do you want ? \n phi =")
          phi=parse(Float64,input())
          Xout=cos.(omega*Xin+phi)
          XT=("cos($omega x+$phi)")
        elseif TrigN=="tan"
          pritnln("What Angulat fequency (omega) do you want ? \n omega =")
          omega=parse(Float64,input())
          println("What Dephasage (phi) do you want ? \n phi =")
          phi=parse(Float64,input())
          Xout=tan.(omega*Xin+phi)
          XT=("tan($omega x+$phi)")
        elseif TrigN=="cot"
          pritnln("What Angulat fequency (omega) do you want ? \n omega =")
          omega=parse(Float64,input())
          println("What Dephasage (phi) do you want ? \n phi =")
          phi=parse(Float64,input())
          Xout=cot.(omega*Xin+phi)
          XT=("cot($omega x+$phi)")
        end
      elseif TTrig=="Inverse"
        println("What type do you want ? (Possible choices: acos, asin, atan or acot) \n Your choice is :")
        TrigN=input()
        if TrigN=="asin"
          pritnln("Vector is standardize (Divided by Max)")
          println("f ? : asin(XStd/f) (abs(f)>1)? \n f=")
          f=parse(Float64,input())
          if abs(f)<1
            f=1
          end
          Xout=asin.(Xin/maximum(abs(Xin))/f)
          XT=("asin(XStd/$f)")
        elseif TrigN=="acos"
          pritnln("Vector is standardize (Divided by Max)")
          println("f ? : acos(XStd/f) (abs(f)>1)? \n f=")
          f=parse(Float64,input())
          if abs(f)<1
            f=1
          end
          Xout=acos.(Xin/maximum(abs(Xin))/f)
          XT=("acos(XStd/$f)")
        elseif TrigN=="atan"
          println("What Factor f ? (atan(fX+t)) \n f=")
          f=parse(Float64,input())
          println("What terme ? \n t=")
          t=parse(Float64,input())
          Xout=atan.(f*Xin+t)
          XT=("atan(fX+t)")
        elseif TrigN=="acot"
          println("What Factor f ? (acot(fX+t)) \n f=")
          f=parse(Float64,input())
          println("What terme ? \n t=")
          t=parse(Float64,input())
          Xout=acot.(f*Xin+t)
          XT=("acot(fX+t)")
        end
      end
    elseif SwCu=="Frac"
      println("Order of the wanted term : y=x^-{Order} ? \n Order=")
      u=parse(Float64,input())
      if isinteger(u)
        Xu=Xin
      else
        Xu=Xin
        Xu[find(Xu.<=0)]=ones(size(find(Xu.<0)))
        println("All smaller than or equal to 0 elements have been replaced by one for Frac Order $u")
      end
      Xout=1./(Xu.^(u))
      XT=("Frac Order $u")
    elseif SwCu=="Expo"
      println("Order of the wanted Exponential ? : y=exp(x*{Order}) ? \n Order=")
      u=parse(Float64,input())
      Xout=exp.(u*Xin)
      XT=("exp($u x)")
    elseif SwCu=="Loga"
      println("Terme of the wanted Logarithme ? : y=log|x+t|) ? \n t=")
      t=parse(Float64,input())
      Xout=log.(abs(Xin+t))
      XT=("log(x+t)")
    end
  end
return Xout,XT
end

###############################################################################

# Statistics discrabing the Regression using ANOVA analysis thanks to distributions
# functions

function StatReg(Xin,Yin,XMat,XName,p)
  Yb=mean(Yin)
  SquareX=(transpose(XMat)*XMat)
  n=length(Yin)
  Bh=(SquareX^-1)*transpose(XMat)*Yin # Estimator for Inv(Mat) ?
  Yh=XMat*Bh
  SSE=sum((Yin-Yh).^2)
  SSR=sum((Yh-Yb).^2)
  SSTot=sum((Yin-Yb).^2)
  R2=SSR/SSTot # the biggest the best
  dfR=p-1
  MSR=SSR/dfR
  dfE=n-p
  MSE=SSE/dfE
  F=MSE/MSR
  # ---------------------- Calc p-Value Model --------------------------
  FD=FDist(dfR,dfE)
  pModel=cdf(FD,F) #Smallest the best
  # p- Values indie
  TD=TDist(dfE)
  pVect=ones(size(Bh))
  for i=1:length(Bh)
    Yh=XMat[:,i]*Bh[i]
    SSE=sum((Yin-Yh).^2)
    if XName[i]=="Constant"
      SX2=sum((Xin-mean(Xin)).^2)
      SEBh=sqrt(SSE*(1/n+mean(Xin)^2/SX2))
    else
      X=XMat[:,i]
      SX2=sum((X-mean(X)).^2)
      SEBh=sqrt(SSE/SX2)
    end
    t=Bh[i]*sqrt(dfE)/SEBh
    pVar=2*(1-cdf(TD,t))
    pVect[i]=pVar
  end

  return pVect,pModel,Bh
end

#############################################################################

# Main function of Multiple Regression implementing ANOVA tests
# and a user friendly approach to regression

function MultiRegOpt(Xin::Vector, Yin::Vector, FunFamily="Poly"::String, p=3::Int64,Stepwise="None"::String, threshold=0.05::Float64)
  n=length(Xin)
  pVModels=[]
  XMatArray=[]
  BArray=[]
  YhA=[]
  XMatA=[]
  pVarA=[]
  XNameA=[]
  NbrModelTested=0
  if FunFamily=="Custom"
    NbrModel=p
    XMat=ones(Float64,n,NbrModel)
    XName=[]
    for Counti=0:NbrModel-1 # Nombre of function involved
      if Counti==0
        FX=ones(size(Xin))
        push!(XName,"Constant")
      else
        Xout=NextFun(Xin,FunFamily,Counti)
        FX=Xout[1]
        push!(XName,Xout[2])
      end
      XMat[:,Counti+1]=FX
    end
    pVect,pModel,Bh=StatReg(Xin,Yin,XMat,XName,p)
    # Stepwise, taking out not necessary variables
    OldXMat=XMat
    OldpVect=pVect
    OldBh=Bh
    OldXMat=XMat
    OldpModel=pModel
    while sum(pVect.>threshold)>0 && Stepwise=="Student"# There is pVect>0.05
      i=findmax(pVect)[2] #eliminated fun
      if (NbrModel-1)==0
        println("Impossible to find a acceptable regression with given number of function for threshold=$threshold, First Result given")
        XMat=OldXMat
        pVect=OldpVect
        Bh=OldBh
        XMat=OldXMat
        break
      else
        NewXMat=ones(Float64,n,NbrModel-1)
        NewXName=[]
        for j=1:NbrModel
          if j<i
            NewXMat[:,j]=XMat[:,j]
            push!(NewXName,XName[j])
          elseif j>i
            NewXMat[:,j-1]=XMat[:,j]
            push!(NewXName,XName[j])
          end
        end
        println("MultiReg New:NewXName=$NewXName")
        pVect,pModel,Bh=StatReg(Xin,Yin,NewXMat,NewXName,p)
        XMat=NewXMat
        XName=NewXName
        NbrModel=NbrModel-1
      end
    end
    push!(pVarA,pVect)
    push!(pVModels,pModel)
    push!(BArray,Bh)
    push!(XMatA,XMat)
    push!(XNameA,XName)
  else
    for NbrModel=1:p
      XMat=ones(Float64,n,NbrModel)
      XName=[]
      for Counti=0:NbrModel-1 # Nombre of function involved
        if Counti==0
          FX=ones(size(Xin))
          push!(XName,"Constant")
        else
          Xout=NextFun(Xin,FunFamily,Counti)
          FX=Xout[1]
          push!(XName,Xout[2])
        end
        XMat[:,Counti+1]=FX
      end
      pVect,pModel,Bh=StatReg(Xin,Yin,XMat,XName,p)
      # Stepwise, taking out not necessary variables
      OldpVect=pVect
      OldBh=Bh
      OldXMat=XMat
      OldpModel=pModel
      while sum(pVect.>threshold)>0 && Stepwise=="Student"# there is pVect>0.05
        i=findmax(pVect)[2] #eliminated fun
        if (NbrModel-1)==0
          println("Impossible to find a acceptable regression with given number of function for threshold=$threshold, First Result given")
          XMat=OldXMat
          pVect=OldpVect
          Bh=OldBh
          XMat=OldXMat
          break
        else
          NewXMat=ones(Float64,n,NbrModel-1)
          NewXName=[]
          for j=1:NbrModel
            if j<i
              NewXMat[:,j]=XMat[:,j]
              push!(NewXName,XName[j])
            elseif j>i
              NewXMat[:,j-1]=XMat[:,j]
              push!(NewXName,XName[j])
            end
          end
          pVect,pModel,Bh=StatReg(Xin,Yin,NewXMat,NewXName,p)
          XMat=NewXMat
          XName=NewXName
          NbrModel=NbrModel-1
        end
      end
      if Stepwise=="Fleussu"
        Yb=mean(Yin)
        if sum(abs(Bh).<Yb*10.0^-5)>0
          println(abs(Bh).<Yb*10.0^-5)
          Bh[abs(Bh).<Yb*10.0^-5]=maximum(Bh)*2
          i=findin(X,maximum(X))
          NewXName=[]
          NewBh=[]
          for j=1:length(Bh)
            if sum(findin(i,j))==0
              push!(NewXName,XName[j])
              push!(NewBh,Bh[j])
            end
          end
          XName=NewXName
          Bh=NewBh
          NbrModel=length(Bh)
        end
      end
      push!(pVarA,pVect)
      push!(pVModels,pModel)
      push!(BArray,Bh)
      push!(XMatA,XMat)
      push!(XNameA,XName)
    end
  end
  return BArray,pVModels,XMatA,pVarA,XNameA
end

############################################################################
end # module


text:
using Documenter
using ConstrainedDynamicsVis

makedocs(;
    modules = [ConstrainedDynamicsVis],
    format = Documenter.HTML(
        canonical = "https://janbruedigam.github.io/ConstrainedDynamicsVis.jl/stable/",
        assets = ["assets/favicon.ico"],
    ),
    pages = [
        "Home" => "index.md",
    ],
    sitename = "ConstrainedDynamicsVis.jl"
)

deploydocs(; repo = "github.com/janbruedigam/ConstrainedDynamicsVis.jl")

ext: jl


text:
# This file was generated, do not modify it.

using MLJ

@load LinearRegressor pkg=MLJLinearModels

using RDatasets, DataFrames
boston = dataset("MASS", "Boston")
first(boston, 3)

describe(boston, :mean, :std, :eltype)

using ScientificTypes
data = coerce(boston, :Tax=>Continuous, :Rad=>Continuous);

y = data.MedV
X = select(data, Not(:MedV));

mdl = LinearRegressor()

mach = machine(mdl, X, y)
fit!(mach)

fp = fitted_params(mach)
@show round.(fp.coefs[1:3], sigdigits=3)
@show round(fp.intercept, sigdigits=3)

ŷ = predict(mach, X)
round(rms(ŷ, y), sigdigits=4)

X2 = hcat(X, X.LStat .* X.Age);

rename!(X2, :x1 => :interaction);

mach = machine(mdl, X2, y)
fit!(mach)
ŷ = predict(mach, X2)
round(rms(ŷ, y), sigdigits=4)

X3 = hcat(X.LStat, X.LStat.^2)
machine(mdl, X3, y)
fit!(mach)
ŷ = predict(mach, X3)
round(rms(ŷ, y), sigdigits=4)

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
#=
cudatest:
- Julia version:
- Author: jerrymei
- Date: 2018-05-29
=#
using CUDAdrv
using CUDAnative
import NVTX
import CUDNN

NVTX.stop() # Use `--profile-from-start=off`
ITERATIONS = 5

img=CuArray(rand(Float32,144,144,64,64))
fil=CuArray(rand(Float32,3,3,64,64))
dest=CuArray(rand(Float32,144,144,64,64))

# CUDAnative.@device_code_typed @cuda (10,10,shmem) fil_trans_kernel(fil,fil_trans)

function benchmark_fil_trans()
    NVTX.@range "kernel" begin
        CUDNN.conv2d!(dest,img,fil)
        CUDAdrv.synchronize()
    end
end

#warmup
NVTX.mark("Warmup benchmark")
benchmark_fil_trans()

NVTX.@activate CUDAdrv.@profile begin
    for i in 1:ITERATIONS
        benchmark_fil_trans()
    end
end

max_stars_repo_name: JerryLingjieMei/WinogradConvolution

max_issues_repo_name: JerryLingjieMei/WinogradConvolution

max_forks_repo_name: JerryLingjieMei/WinogradConvolution


text:
abstract type TabuStrategy end

struct NoTabu <: TabuStrategy end

struct KeenTabu <: TabuStrategy
    tabu_tenure::Int
    tabu_list::Dictionary{Int,Int}
end

struct WeakTabu <: TabuStrategy
    tabu_tenure::Int
    pick_tenure::Int
    tabu_list::Dictionary{Int,Int}
end

tabu() = NoTabu()
function tabu(tabu_tenure)
    tabu_list = Dictionary{Int, Int}()
    return KeenTabu(tabu_tenure, tabu_list)
end
function tabu(tabu_tenure, pick_tenure)
    tabu_list = Dictionary{Int, Int}()
    return WeakTabu(tabu_tenure, pick_tenure, tabu_list)
end

tenure(strategy, ::Val{:tabu}) = strategy.tabu_tenure
tenure(strategy::WeakTabu, ::Val{:pick}) = strategy.pick_tenure
tenure(::TabuStrategy, ::Val{:pick}) = zero(Int)
tenure(::NoTabu, ::Val{:tabu}) = zero(Int)
tenure(strategy, field) = tenure(strategy, Val(field))

"""
    _tabu(s::S) where S <: Union{_State, AbstractSolver}
Access the list of tabu variables.
"""
tabu_list(ts) = ts.tabu_list

"""
    _tabu(s::S, x) where S <: Union{_State, AbstractSolver}
Return the tabu value of variable `x`.
"""
tabu_value(ts, x) = tabu_list(ts)[x]

"""
    _decrease_tabu!(s::S, x) where S <: Union{_State, AbstractSolver}
Decrement the tabu value of variable `x`.
"""
decrease_tabu!(ts, x) = tabu_list(ts)[x] -= 1

"""
    _delete_tabu!(s::S, x) where S <: Union{_State, AbstractSolver}
Delete the tabu entry of variable `x`.
"""
delete_tabu!(ts, x) = delete!(tabu_list(ts), x)

"""
    _empty_tabu!(s::S) where S <: Union{_State, AbstractSolver}
Empty the tabu list.
"""
empty_tabu!(ts) = Dictionaries.empty!(tabu_list(ts))

"""
    _length_tabu!(s::S) where S <: Union{_State, AbstractSolver}
Return the length of the tabu list.
"""
length_tabu(ts) = length(tabu_list(ts))

"""
    _insert_tabu!(s::S, x, tabu_time) where S <: Union{_State, AbstractSolver}
Insert the bariable `x` as tabu for `tabu_time`.
"""
insert_tabu!(ts, x, ::Val{:tabu}) = insert!(tabu_list(ts), x, max(1, tenure(ts, :tabu)))
insert_tabu!(ts::KeenTabu, x, kind::Symbol) = insert_tabu!(ts, x, Val(kind))
insert_tabu!(ts::WeakTabu, x, kind) = insert!(tabu_list(ts), x, max(1, tenure(ts, kind)))
insert_tabu!(ts, x, kind) = nothing

"""
    _decay_tabu!(s::S) where S <: Union{_State, AbstractSolver}
Decay the tabu list.
"""
function decay_tabu!(ts)
    foreach(
        ((x, tabu),) -> tabu == 1 ? delete_tabu!(ts, x) : decrease_tabu!(ts, x),
        pairs(tabu_list(ts))
    )
end


text:
export disp, showdict, log, @log, setfilelogging, setlogfile, errorstring

disp(x...) = println(join(map(string,x),", "))

function showdict(a, indent = "")
    sk = sort(collect(keys(a)))
    for k in sk
        print("$(indent)$k: ")
        if isa(a[k], Dict)
            println("")
            showdict(a[k], indent*"  ")
        else
            println("$(a[k])")
        end
    end
end

function setfilelogging(a::Bool)
    global LOGTOFILE
    LOGTOFILE = a
end

function setlogfile(a::AbstractString)
    global LOGFILE
    LOGFILE = a
end

import Base.log
log(a::AbstractString, args...; kargs...) = log(STDOUT, a, args...; kargs...)
function log(io::IO, a::AbstractString, args...; indent = 0, tofile = [], toSTDOUT = true)
    buf = IOBuffer()
    println(buf, Libc.strftime("%Y-%m-%d %T %z %Z", time()), "  |  ", repeat("  ", indent), join([a, args...], " "))
    str = takebuf_string(buf)

    toSTDOUT && println(io, str[end] == '\n' ? str[1:end-1] : str)
    if (tofile == [] && LOGTOFILE) || tofile == true || isa(tofile, AbstractString)
        filename = isa(tofile, AbstractString) ? tofile : LOGFILE
        @spawnat 1 open(filename, "a") do fid
            write(fid, str)
        end
    end
    nothing
end

macro log(a...)
    quote
        log($(a...))
    end
end

function errorstring(e)
    buf = IOBuffer()
    showerror(buf, e, catch_backtrace())
    takebuf_string(buf)
end


text:
"""
    subs

Substitute values into a symbolic expression.

Examples
```
@vars x y
ex = x^2 + y^2
subs(ex, x, 1) # 1 + y^2
subs(ex, (x, 1)) # 1 + y^2
subs(ex, (x, 1), (y,x)) # 1 + x^2, values are substituted left to right.
subs(ex, x=>1)  # alternate to subs(x, (x,1))
subs(ex, x=>1, y=>1) # ditto
```
"""
function subs(ex::T, var::S, val) where {T<:SymbolicType, S<:SymbolicType}
    s = Basic()
    ccall((:basic_subs2, libsymengine), Nothing, (Ref{Basic}, Ref{Basic}, Ref{Basic}, Ref{Basic}), s, ex, var, val)
    return s
end
function subs(ex::T, d::CMapBasicBasic) where T<:SymbolicType
    s = Basic()
    ccall((:basic_subs, libsymengine), Nothing, (Ref{Basic}, Ref{Basic}, Ptr{Cvoid}), s, ex, d.ptr)
    return s
end

subs(ex::T, d::Dict) where {T<:SymbolicType} = subs(ex, CMapBasicBasic(d))
subs(ex::T, y::Tuple{S, Any}) where {T <: SymbolicType, S<:SymbolicType} = subs(ex, y[1], y[2])
subs(ex::T, y::Tuple{S, Any}, args...) where {T <: SymbolicType, S<:SymbolicType} = subs(subs(ex, y), args...)
subs(ex::T, d::Pair...) where {T <: SymbolicType} = subs(ex, [(p.first, p.second) for p in d]...)

## Allow an expression to be called, as with ex(2). When there is more than one symbol, one can rely on order of `free_symbols` or
## be explicit by passing in pairs : `ex(x=>1, y=>2)` or a dict `ex(Dict(x=>1, y=>2))`.
function (ex::Basic)(args...)
  xs = free_symbols(ex)
  subs(ex, collect(zip(xs, args))...)
end
(ex::Basic)(x::Dict) = subs(ex, x)
(ex::Basic)(x::Pair...) = subs(ex, x...)

## Lambdify

## Map symengine classes to function names
fn_map = Dict(
              :Add => :+,
              :Sub => :-,
              :Mul => :*,
              :Div => :/,
              :Pow => :^,
              :re  => :real,
              :im  => :imag,
              :Abs => :abs # not really needed as default in map_fn covers this
              )

map_fn(key, fn_map) = haskey(fn_map, key) ? fn_map[key] : Symbol(lowercase(string(key)))

function convert(::Type{Expr}, ex::Basic)
    fn = get_symengine_class(ex)

    if fn == :Symbol
        return Symbol(toString(ex))
    elseif (fn in number_types) || (fn == :Constant)
        return N(ex)
    end

    as = get_args(ex)

    Expr(:call, map_fn(fn, fn_map), [convert(Expr,a) for a in as]...)
end

function convert(::Type{Expr}, m::AbstractArray{Basic, 2})
    col_args = []
    for i = 1:size(m, 1)
        row_args = []
        for j = 1:size(m, 2)
            push!(row_args, convert(Expr, m[i, j]))
        end
        row = Expr(:hcat, row_args...)
        push!(col_args, row)
    end
    Expr(:vcat, col_args...)
end

function convert(::Type{Expr}, m::AbstractArray{Basic, 1})
    row_args = []
    for j = 1:size(m, 1)
        push!(row_args, convert(Expr, m[j]))
    end
    Expr(:vcat, row_args...)
end

walk_expression(b) = convert(Expr, b)

"""
    lambdify
evaluates a symbolless expression or returns a function
"""
function lambdify(ex, vars=[])
    if length(vars) == 0
        vars = free_symbols(ex)
    end
    body = convert(Expr, ex)
    lambdify(body, vars)
end

lambdify(ex::BasicType, vars=[]) = lambdify(Basic(ex), vars)

function lambdify(ex::Expr, vars)
    if length(vars) == 0
        # return a number
        eval(ex)
    else
        # return a function
        _lambdify(ex, vars)
    end
end

function _lambdify(ex::Expr, vars)
    try
        fn = eval(Expr(:function,
                  Expr(:call, gensym(), map(Symbol,vars)...),
                       ex))
        (args...) -> invokelatest(fn, args...) # https://github.com/JuliaLang/julia/pull/19784
    catch err
        throw(ArgumentError("Expression does not lambdify"))
    end
end


text:
"""
# Description
Takes NzxNt arrays of profiles for variables u, v, T and returns
Nzx(Nt-1) arrays of the profile evolution for u, v, T, u'w', v'w', and w'T'
the horizontally averaged flux for variable V.

# Arguments
Unscaled u, v, T, z, t, and f
"""
function reconstruct_flux_profiles(u, v, T, z, t, f)

    Δz = diff(z)
    Δt = diff(t, dims=1)'

    Nz,Nt = size(T)

    dudt = (u[:,2:Nt] .- u[:,1:Nt-1]) ./ Δt # Nz x (Nt-1) array of approximate dVdt values
    dvdt = (v[:,2:Nt] .- v[:,1:Nt-1]) ./ Δt # Nz x (Nt-1) array of approximate dVdt values
    dTdt = (T[:,2:Nt] .- T[:,1:Nt-1]) ./ Δt # Nz x (Nt-1) array of approximate dVdt values
    u = u[:,1:Nt-1]
    v = v[:,1:Nt-1]
    T = T[:,1:Nt-1]

    """ evaluates wϕ = ∫ ∂z(wϕ) dz """
    function wϕ(∂z_wϕ)
        ans = zeros(Nz+1, Nt-1) # one fewer column than T
        for i in 1:Nt-1, h in 1:Nz-1
            ans[h+1, i] = ans[h, i] + Δz[h] * ∂z_wϕ[h, i]
        end
        return ans
    end

    duw_dz = -dudt .+ f*v
    dvw_dz = -dvdt .- f*u
    dwT_dz = -dTdt

    # println(size(wV(duw_dz)))
    # u, v, T, uw, vw, wT, t
    return (u, v, T, wϕ(duw_dz), wϕ(dvw_dz), wϕ(dwT_dz), t[1:Nt-1])
end

using ClimateParameterizations, Plots
#
𝒟 = ClimateParameterizations.Data.data("strong_wind")

𝒟_reconstructed = ClimateParameterizations.Data.data("strong_wind_weak_heating", reconstruct_fluxes=true)

𝒟_reconstructed
z = 𝒟_reconstructed.uw.z
t = 𝒟_reconstructed.t
Nt = length(𝒟.t)
output_gif_directory = "TestReconstructFluxes"
animate_gif((𝒟_reconstructed.uw.coarse, 𝒟.uw.coarse[:,1:Nt-1]), z, t, "uw",
            x_label=["∫(-du/dt + fv)dz", "truth"],
            filename="uw_reconstructed",
            directory=output_gif_directory)
animate_gif((𝒟_reconstructed.vw.coarse, 𝒟.vw.coarse[:,1:Nt-1]), z, t, "vw",
            x_label=["∫(-dv/dt - fu)dz", "truth"],
            filename="vw_reconstructed",
            directory=output_gif_directory)
animate_gif((𝒟_reconstructed.wT.coarse, 𝒟.wT.coarse[:,1:Nt-1]), z, t, "wT",
            x_label=["∫(-dw/dt)dz", "truth"],
            filename="wT_reconstructed",
            directory=output_gif_directory)


text:
using Oceananigans.Operators

"""
    BuoyancyTracer <: AbstractBuoyancy{Nothing}

Type indicating that the tracer `b` represents buoyancy.
"""
struct BuoyancyTracer <: AbstractBuoyancy{Nothing} end

required_tracers(::BuoyancyTracer) = (:b,)

@inline buoyancy_perturbation(i, j, k, grid, ::BuoyancyTracer, C) = @inbounds C.b[i, j, k]

@inline ∂x_b(i, j, k, grid, ::BuoyancyTracer, C) = ∂xᶠᵃᵃ(i, j, k, grid, C.b)
@inline ∂y_b(i, j, k, grid, ::BuoyancyTracer, C) = ∂yᵃᶠᵃ(i, j, k, grid, C.b)
@inline ∂z_b(i, j, k, grid, ::BuoyancyTracer, C) = ∂zᵃᵃᶠ(i, j, k, grid, C.b)


text:
function run(s)
    total = 0
    for line in split(s, "\n")
        mass = parse(Int, line)
        fuelneeded = 0
        while fuelneeded >= 0
            total += fuelneeded
            fuelneeded = mass ÷ 3 - 2
            mass = fuelneeded
        end
    end
    return total
end

function main()
    res, time, memory = @timed run(ARGS[1])
    # println(time)
    println("_duration:$(time*1000)")
    println(res)
end

main()


text:
const USE_GPU = false  # Use GPU? If this is set false, then no GPU needs to be available
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 2)
else
    @init_parallel_stencil(Threads, Float64, 2)
end
using Plots, Printf, Statistics, LinearAlgebra

@parallel function compute_timesteps!(dτVx::Data.Array, dτVy::Data.Array, dτPt::Data.Array, Mus::Data.Array, Vsc::Data.Number, Ptsc::Data.Number, min_dxy2::Data.Number, max_nxy::Int)
    @all(dτVx) = Vsc*min_dxy2/@av_xi(Mus)/4.1
    @all(dτVy) = Vsc*min_dxy2/@av_yi(Mus)/4.1
    @all(dτPt) = Ptsc*4.1*@all(Mus)/max_nxy
    return
end

@parallel function compute_P!(∇V::Data.Array, Pt::Data.Array, Vx::Data.Array, Vy::Data.Array, dτPt::Data.Array, dx::Data.Number, dy::Data.Number)
    @all(∇V)  = @d_xa(Vx)/dx + @d_ya(Vy)/dy
    @all(Pt)  = @all(Pt) - @all(dτPt)*@all(∇V)
    return
end

@parallel function compute_τ!(∇V::Data.Array, τxx::Data.Array, τyy::Data.Array, τxy::Data.Array, Vx::Data.Array, Vy::Data.Array,  Mus::Data.Array, dx::Data.Number, dy::Data.Number)
    @all(τxx) = 2.0*@all(Mus)*(@d_xa(Vx)/dx - 1.0/3.0*@all(∇V))
    @all(τyy) = 2.0*@all(Mus)*(@d_ya(Vy)/dy - 1.0/3.0*@all(∇V))
    @all(τxy) = 2.0*@av(Mus)*(0.5*(@d_yi(Vx)/dy + @d_xi(Vy)/dx))
    return
end

@parallel function compute_dV!(Rx::Data.Array, Ry::Data.Array, dVxdτ::Data.Array, dVydτ::Data.Array, Pt::Data.Array, Rog::Data.Array, τxx::Data.Array, τyy::Data.Array, τxy::Data.Array, dampX::Data.Number, dampY::Data.Number, dx::Data.Number, dy::Data.Number)
    @all(Rx)    = @d_xi(τxx)/dx + @d_ya(τxy)/dy - @d_xi(Pt)/dx
    @all(Ry)    = @d_yi(τyy)/dy + @d_xa(τxy)/dx - @d_yi(Pt)/dy + @av_yi(Rog)
    @all(dVxdτ) = dampX*@all(dVxdτ) + @all(Rx)
    @all(dVydτ) = dampY*@all(dVydτ) + @all(Ry)
    return
end

@parallel function compute_V!(Vx::Data.Array, Vy::Data.Array, dVxdτ::Data.Array, dVydτ::Data.Array, dτVx::Data.Array, dτVy::Data.Array)
    @inn(Vx) = @inn(Vx) + @all(dτVx)*@all(dVxdτ)
    @inn(Vy) = @inn(Vy) + @all(dτVy)*@all(dVydτ)
    return
end

@parallel_indices (ix,iy) function bc_x!(A::Data.Array)
    A[1  , iy] = A[2    , iy]
    A[end, iy] = A[end-1, iy]
    return
end

@parallel_indices (ix,iy) function bc_y!(A::Data.Array)
    A[ix, 1  ] = A[ix, 2    ]
    A[ix, end] = A[ix, end-1]
    return
end

##################################################
@views function Stokes2D()
    # Physics
    lx, ly    = 10.0, 10.0  # domain extends
    μs0       = 1.0         # matrix viscosity
    μsi       = 0.1         # inclusion viscosity
    ρgi       = 1.0         # inclusion density*gravity perturbation
    # Numerics
    iterMax   = 10000       # maximum number of pseudo-transient iterations
    nout      = 200         # error checking frequency
    Vdmp      = 4.0         # damping paramter for the momentum equations
    Vsc       = 1.0         # relaxation paramter for the momentum equations pseudo-timesteps limiters
    Ptsc      = 1.0/4.0     # relaxation paramter for the pressure equation pseudo-timestep limiter
    ε         = 1e-6        # nonlinear absolute tolerence
    nx, ny    = 127, 127    # numerical grid resolution; should be a mulitple of 32-1 for optimal GPU perf
    # Derived numerics
    dx, dy    = lx/(nx-1), ly/(ny-1) # cell sizes
    min_dxy2  = min(dx,dy)^2
    max_nxy   = max(nx,ny)
    dampX     = 1.0-Vdmp/nx # damping term for the x-momentum equation
    dampY     = 1.0-Vdmp/ny # damping term for the y-momentum equation
    # Array allocations
    Pt        = @zeros(nx  ,ny  )
    dτPt      = @zeros(nx  ,ny  )
    ∇V        = @zeros(nx  ,ny  )
    Vx        = @zeros(nx+1,ny  )
    Vy        = @zeros(nx  ,ny+1)
    τxx       = @zeros(nx  ,ny  )
    τyy       = @zeros(nx  ,ny  )
    τxy       = @zeros(nx-1,ny-1)
    Rx        = @zeros(nx-1,ny-2)
    Ry        = @zeros(nx-2,ny-1)
    dVxdτ     = @zeros(nx-1,ny-2)
    dVydτ     = @zeros(nx-2,ny-1)
    dτVx      = @zeros(nx-1,ny-2)
    dτVy      = @zeros(nx-2,ny-1)
    # Initial conditions
    Radc      =  zeros(nx  ,ny  )
    Rog       =  zeros(nx  ,ny  )
    Mus       = μs0*ones(nx,ny)
    Radc     .= [((ix-1)*dx-0.5*lx)^2 + ((iy-1)*dy-0.5*ly)^2 for ix=1:size(Radc,1), iy=1:size(Radc,2)]
    Mus[Radc.<1.0] .= μsi
    Rog[Radc.<1.0] .= ρgi
    Mus       = Data.Array(Mus)
    Rog       = Data.Array(Rog)
    # Preparation of visualisation
    ENV["GKSwstype"]="nul"; if isdir("viz2D_out")==false mkdir("viz2D_out") end; loadpath = "./viz2D_out/"; anim = Animation(loadpath,String[])
    println("Animation directory: $(anim.dir)")
    X, Y, Yv  = 0:dx:lx, 0:dy:ly, (-dy/2):dy:(ly+dy/2)
    # Time loop
    @parallel compute_timesteps!(dτVx, dτVy, dτPt, Mus, Vsc, Ptsc, min_dxy2, max_nxy)
    err=2*ε; iter=1; niter=0; err_evo1=[]; err_evo2=[]
    while err > ε && iter <= iterMax
        if (iter==11)  global wtime0 = Base.time()  end
        @parallel compute_P!(∇V, Pt, Vx, Vy, dτPt, dx, dy)
        @parallel compute_τ!(∇V, τxx, τyy, τxy, Vx, Vy, Mus, dx, dy)
        @parallel compute_dV!(Rx, Ry, dVxdτ, dVydτ, Pt, Rog, τxx, τyy, τxy, dampX, dampY, dx, dy)
        @parallel compute_V!(Vx, Vy, dVxdτ, dVydτ, dτVx, dτVy)
        @parallel (1:size(Vx,1), 1:size(Vx,2)) bc_y!(Vx)
        @parallel (1:size(Vy,1), 1:size(Vy,2)) bc_x!(Vy)
        if mod(iter,nout)==0
            global mean_Rx, mean_Ry, mean_∇V
            mean_Rx = mean(abs.(Rx)); mean_Ry = mean(abs.(Ry)); mean_∇V = mean(abs.(∇V))
            err = maximum([mean_Rx, mean_Ry, mean_∇V])
            push!(err_evo1, maximum([mean_Rx, mean_Ry, mean_∇V])); push!(err_evo2,iter)
            @printf("Total steps = %d, err = %1.3e [mean_Rx=%1.3e, mean_Ry=%1.3e, mean_∇V=%1.3e] \n", iter, err, mean_Rx, mean_Ry, mean_∇V)
        end
        iter+=1; niter+=1
    end
    # Performance
    wtime    = Base.time() - wtime0
    A_eff    = (3*2)/1e9*nx*ny*sizeof(Data.Number)  # Effective main memory access per iteration [GB] (Lower bound of required memory access: Te has to be read and written: 2 whole-array memaccess; Ci has to be read: : 1 whole-array memaccess)
    wtime_it = wtime/(niter-10)                     # Execution time per iteration [s]
    T_eff    = A_eff/wtime_it                       # Effective memory throughput [GB/s]
    @printf("Total steps = %d, err = %1.3e, time = %1.3e sec (@ T_eff = %1.2f GB/s) \n", niter, err, wtime, round(T_eff, sigdigits=2))
    # Visualisation
    p1 = heatmap(X,  Y, Array(Pt)', aspect_ratio=1, xlims=(X[1],X[end]), ylims=(Y[1],Y[end]), c=:inferno, title="Pressure")
    p2 = heatmap(X, Yv, Array(Vy)', aspect_ratio=1, xlims=(X[1],X[end]), ylims=(Yv[1],Yv[end]), c=:inferno, title="Vy")
    p4 = heatmap(X[2:end-1], Yv[2:end-1], log10.(abs.(Array(Ry)')), aspect_ratio=1, xlims=(X[2],X[end-1]), ylims=(Yv[2],Yv[end-1]), c=:inferno, title="log10(Ry)")
    p5 = plot(err_evo2,err_evo1, legend=false, xlabel="# iterations", ylabel="log10(error)", linewidth=2, markershape=:circle, markersize=3, labels="max(error)", yaxis=:log10)
    # display(plot(p1, p2, p4, p5))
    plot(p1, p2, p4, p5); frame(anim)
    gif(anim, "Stokes2D.gif", fps = 15)
    return
end

Stokes2D()


text:
module ConfParser

import Base: haskey, merge!

export ConfParse, parse_conf!, erase!, save!, retrieve, commit!, haskey, merge!

Base.@deprecate open_fh(filename::String, mode::String) open(filename, mode) false

mutable struct ConfParse
    _fh::IO
    _filename::String
    _syntax::String
    _data::Dict
    _is_modified::Bool

    function ConfParse(filename::AbstractString, syntax::AbstractString="")
        isempty(filename) && throw(ArgumentError("no file name specified"))
        fh = open(filename, "r")
        if isempty(syntax)
            syntax = guess_syntax(fh)
        elseif !(syntax == "ini" || syntax == "http" || syntax == "simple" || syntax == "keyonly")
            close(fh)
            throw(ArgumentError("unrecognized configuration syntax: $syntax"))
        end
        conf = new(fh, filename, syntax, Dict(), false)
        finalizer(c->close(c._fh), conf)
        return conf
    end
end

"""
    guess_syntax(fh::IO)

Attempts to guess the configuration file syntax using
regular expressions.
"""
function guess_syntax(fh::IO)
    cnt = countlines(fh)
    seekstart(fh)
    if cnt == 1
        value = String(strip(read(fh, String)))
        all(isxdigit, value) && return "keyonly"
    else
        for line in eachline(fh)
            # is a commented line
            occursin(r"^\s*(?:#|$)", line) && continue

            # is not alphanumeric
            occursin(r"\w", line) || continue

            # remove \n
            line = chomp(line)

            # contains a [block]; ini
            occursin(r"^\s*\[\s*[^\]]+\s*\]\s*$", line) && return "ini"

            # key/value pairs are seperated by a '='; ini
            occursin(r"^\s*[\w-]+\s*=\s*.*\s*$", line) && return "ini"

            # key/value pairs are seperated by a ':'; http
            occursin(r"^\s*[\w-]+\s*:\s*.*\s*$", line) && return "http"

            # key/value pairs are seperated by whitespace; simple
            occursin(r"^\s*[\w-]+\s+.*$", line) && return "simple"
        end
    end

    error("unable to identify the configuration file syntax")
end

"""
    parse_conf!(s::ConfParse)

Tasks appropriate parser function based on configuration syntax.
"""
function parse_conf!(s::ConfParse)
    if s._syntax == "ini"
        parse_ini(s)
    elseif s._syntax == "http"
        parse_http(s)
    elseif s._syntax == "simple"
        parse_simple(s)
    elseif s._syntax == "keyonly"
        parse_keyonly(s)
    else
        error("unknown configuration syntax: $(s._syntax)")
    end
end

"""
    parse_line(line::String)

Sperates by commas, removes newlines and such.
"""
function parse_line(line::String)
    parsed   = String[]
    splitted = split(line, ",")
    for raw = splitted
        if occursin(r"\S+", raw)
            clean = match(r"\S+", raw)
            push!(parsed, clean.match)
        end
    end
    parsed
end

"""
    parse_ini(s::ConfParse)

Parses configuration files utilizing ini sytnax.
Populate the ConfParser.data dictionary.
"""
function parse_ini(s::ConfParse)
    blockname = "default"
    seekstart(s._fh)
    for line in eachline(s._fh)
        # skip comments and newlines
        occursin(r"^\s*(\n|\#|;)", line) && continue

        occursin(r"\w", line) || continue

        line = chomp(line)

        # parse blockname
        m = match(r"^\s*\[\s*([^\]]+)\s*\]$", line)
        if m !== nothing
            blockname = lowercase(m.captures[1])
            continue
        end

        # parse key/value
        m = match(r"^\s*([^=]*[^\s])\s*=\s*(.*)\s*$", line)
        if m !== nothing
            key::String, values::String = m.captures
            if !haskey(s._data, blockname)
                s._data[blockname] = Dict(key => parse_line(values))
            else
                merge!(s._data[blockname], Dict(key => parse_line(values)))
            end
            continue
        end

        error("invalid syntax on line: $(line)")
    end
    nothing
end

"""
    parse_http(s::ConfParse)

Parses configuration files utilizing http sytnax.
Populate the ConfParser.data dictionary.
"""
function parse_http(s::ConfParse)
    seekstart(s._fh)
    for line in eachline(s._fh)
        # skip comments and newlines
        occursin(r"^\s*(\n|\#|;)", line) && continue

        occursin(r"\w", line) || continue

        line = chomp(line)

        m = match(r"^\s*([\w-]+)\s*:\s*(.*)$", line)
        if m !== nothing
            key::String, values::String = m.captures
            s._data[key] = parse_line(values)
            continue
        end

        error("invalid syntax on line: $(line)")
    end
end

"""
    parse_simple(s::ConfParse)

Parses configuration files utilizing simple syntax.
Populates the ConfParser.data dictionary.
"""
function parse_simple(s::ConfParse)
    seekstart(s._fh)
    for line in eachline(s._fh)
        # skip comments and newlines
        occursin(r"^\s*(\n|\#|;)", line) && continue

        occursin(r"\w", line) || continue

        line = chomp(line)

        m = match(r"^\s*([\w-]+)\s+(.*)\s*$", line)
        if m !== nothing
            key::String, values::String = m.captures
            s._data[key] = parse_line(values)
            continue
        end

        error("invalid syntax on line: $(line)")
    end
    nothing
end

"""
    parse_keyonly(s::ConfParse)

Parses configuration files utilizing keyonly syntax.
Populates the ConfParser.data dictionary.
"""
function parse_keyonly(s::ConfParse)
    seekstart(s._fh)
    countlines(s._fh) != 1 && error("keyonly syntax requires a single line.")
    seekstart(s._fh)
    value = String(strip(read(s._fh, String)))
    if !all(isxdigit, value)
        error("Invalid characters in keyonly conf: $(value)")
    end
    s._data["key"] = value

    nothing
end

"""
    craft_content(s::ConfParse)

Craft content strings from data array for saved config.
"""
function craft_content(s::ConfParse)
   content = IOBuffer()
   if s._syntax == "ini"
        for (block, key_values) = s._data
            println(content, "[$block]")
            for (key, values) = key_values
                if values isa Vector{<:AbstractString}
                    print(content, key, '=')
                    join(content, values, ',')
                    println(content)
                else
                    println(content, key, '=', values)
                end
            end
            println(content)
        end
    elseif s._syntax == "http"
        for (key, values) = s._data
            if values isa Vector{<:AbstractString}
                print(content, key, ": ")
                join(content, values, ',')
                println(content)
            else
                println(content, key, ": ", values)
            end
        end
    elseif s._syntax == "simple"
        for (key, values) = s._data
            if values isa Vector{<:AbstractString}
                print(content, key, ' ')
                join(content, values, ',')
                println(content)
            else
                println(content, key, ' ', values)
            end
        end
    elseif s._syntax == "keyonly"
        println(content,s._data["key"])
    else
        error("unknown syntax type: $(s._syntax)")
    end
    String(take!(content))
end

"""
    erase!(s::ConfParse, block::String, key::String)

Remove entry from ini block.
"""
function erase!(s::ConfParse, block::String, key::String)
    block_key = getkey(s._data, lowercase(block), nothing)
    if block_key !== nothing
        if haskey(s._data[block_key], key)
            delete!(s._data[block_key], key)
            s._is_modified = true
        end
    end
    s._is_modified
end

"""
    erase!(s::ConfParse, key::String)

Remove entry from config (outside of block if ini).
"""
function erase!(s::ConfParse, key::String)
    if haskey(s._data, key)
        delete!(s._data, key)
        s._is_modified = true
    end
    s._is_modified
end

"""
    save!(s::ConfParse, filename=nothing)

Write out new or modified configuration files.
"""
function save!(s::ConfParse, filename=nothing)
    # if data has not been modified and a new file has not
    # been specified, don't write out
    !s._is_modified && filename === nothing && return

    # if there is no content to write out, don't create an
    # empty file
    content = craft_content(s)
    content === nothing && return

    if filename === nothing
        s._fh = open(s._filename, "w")
    else
        s._fh = open(filename, "w")
    end

    write(s._fh, content)
    flush(s._fh)
    nothing
end

"""
    retrieve(s::ConfParse, key::String)

Retrieve data outside of a block.
"""
function retrieve(s::ConfParse, key::String)
    k = s._data[key]
    length(k) == 1 ? first(k) : k
end

"""
    retrieve(s::ConfParse, block::String, key::String)

Retrieve data from an ini config file block.
"""
function retrieve(s::ConfParse, block::String, key::String)
    k = s._data[lowercase(block)][key]
    length(k) == 1 ? first(k) : k
end

"""
    retrieve(s::ConfParse, key::String, t::Type)

Retrieve data outside of a block and converting to type.
"""
function retrieve(s::ConfParse, key::String, t::Type)
    k = s._data[key]
    parse(t, length(k) == 1 ? first(k) : k)
end

"""
    retrieve(s::ConfParse, block::String, key::String, t::Type)

Retrieve data from an ini config file block and converting to type.
"""
function retrieve(s::ConfParse, block::String, key::String, t::Type)
    k = s._data[lowercase(block)][key]
    parse(t, length(k) == 1 ? first(k) : k)
end

"""
    commit!(s::ConfParse, key::String, value::Any)

Insert data in a config file.
"""
function commit!(s::ConfParse, key::String, value::Any)
    s._data[key]   = value
    s._is_modified = true
end

"""
    commit!(s::ConfParse, block::String, key::String, values::String)

Insert data inside an ini file block.
"""
function commit!(s::ConfParse, block::String, key::String, values::String)
    if s._syntax != "ini"
        throw(ArgumentError("invalid setter function called for syntax type: $(s._syntax)"))
    end
    s._data[lowercase(block)][key] = [values]
    s._is_modified                 = true
end

"""
    Base.merge!(s::ConfParse, t::ConfParse)

Merge data of two configuration files.
"""
function merge!(s::ConfParse, t::ConfParse)
    for (block, key_values) = t._data
        if !haskey(s._data, block)
            s._data[block] = key_values
        else
            merge!(s._data[block], key_values)
        end
    end
    s._is_modified = true
end

"""
    haskey(s::ConfParse, key::String)

Check if a key exists.
"""
haskey(s::ConfParse, key::String) = haskey(s._data, key)

"""
    haskey(s::ConfParse, block::String, key::String)

Check if a key exists inside an ini file block.
"""
haskey(s::ConfParse, block::String, key::String) = haskey(s._data, lowercase(block)) && haskey(s._data[lowercase(block)], key)

end # module ConfParser


text:
module API

import Pkg3
using Pkg3.Types
using Base.Random.UUID

previewmode_info() = info("In preview mode")

add(pkg::String; kwargs...)               = add([pkg]; kwargs...)
add(pkgs::Vector{String}; kwargs...)      = add([PackageSpec(pkg) for pkg in pkgs]; kwargs...)
add(pkgs::Vector{PackageSpec}; kwargs...) = add(EnvCache(), pkgs; kwargs...)

function add(env::EnvCache, pkgs::Vector{PackageSpec}; preview::Bool=env.preview[])
    env.preview[] = preview
    preview && previewmode_info()
    project_resolve!(env, pkgs)
    registry_resolve!(env, pkgs)
    ensure_resolved(env, pkgs, true)
    Pkg3.Operations.add(env, pkgs)
end

rm(pkg::String; kwargs...)               = rm([pkg]; kwargs...)
rm(pkgs::Vector{String}; kwargs...)      = rm([PackageSpec(pkg) for pkg in pkgs]; kwargs...)
rm(pkgs::Vector{PackageSpec}; kwargs...) = rm(EnvCache(), pkgs; kwargs...)

function rm(env::EnvCache, pkgs::Vector{PackageSpec}; preview::Bool=env.preview[])
    env.preview[] = preview
    preview && previewmode_info()
    project_resolve!(env, pkgs)
    manifest_resolve!(env, pkgs)
    Pkg3.Operations.rm(env, pkgs)
end

up(;kwargs...)                           = up(PackageSpec[], kwargs...)
up(pkg::String; kwargs...)               = up([pkg]; kwargs...)
up(pkgs::Vector{String}; kwargs...)      = up([PackageSpec(pkg) for pkg in pkgs]; kwargs...)
up(pkgs::Vector{PackageSpec}; kwargs...) = up(EnvCache(), pkgs; kwargs...)

function up(env::EnvCache, pkgs::Vector{PackageSpec};
            level::UpgradeLevel=UpgradeLevel(:major), mode::Symbol=:project, preview::Bool=env.preview[])
    env.preview[] = preview
    preview && previewmode_info()
    if isempty(pkgs)
        if mode == :project
            for (name::String, uuid::UUID) in env.project["deps"]
                push!(pkgs, PackageSpec(name, uuid, level))
            end
        elseif mode == :manifest
            for (name, infos) in env.manifest, info in infos
                uuid = UUID(info["uuid"])
                push!(pkgs, PackageSpec(name, uuid, level))
            end
        end
    else
        project_resolve!(env, pkgs)
        manifest_resolve!(env, pkgs)
        ensure_resolved(env, pkgs)
    end
    Pkg3.Operations.up(env, pkgs)
end

test(;kwargs...)                           = test(PackageSpec[], kwargs...)
test(pkg::String; kwargs...)               = test([pkg]; kwargs...)
test(pkgs::Vector{String}; kwargs...)      = test([PackageSpec(pkg) for pkg in pkgs]; kwargs...)
test(pkgs::Vector{PackageSpec}; kwargs...) = test(EnvCache(), pkgs; kwargs...)

function test(env::EnvCache, pkgs::Vector{PackageSpec}; coverage=false, preview=env.preview[])
    env.preview[] = preview
    preview && previewmode_info()
    project_resolve!(env, pkgs)
    manifest_resolve!(env, pkgs)
    ensure_resolved(env, pkgs)
    Pkg3.Operations.test(env, pkgs; coverage=coverage)
end

end # module


text:
using MAT
# using PyCall

function generate_MAT_file(rnx::Array, species::Array, stoichiometry::Array)
  filename = "model.mat"
  file = matopen(filename, "w")
  rxns = [rr.rnxName for rr in rnx]
  write(file, "rxns", rxns)
  write(file, "S", stoichiometry)
  lb = zeros(length(rxns), 1)
  write(file, "lb", lb)
  ub = 1000 * ones(length(rxns), 1)
  write(file, "ub", ub)
  c = zeros(length(rxns), 1)
  write(file, "c", c)
  write(file, "mets", species)
  b = 100 * ones(length(species), 1)
  write(file, "b", b)
  write(file, "osenseStr", "max")
  csense = ["L" for i = 1:length(species)]
  write(file, "csense", csense)

  # write(file, "genes", [""])
  # write(file, "rules", ["" for i = 1:length(rxns)])
  close(file)
  return filename
end

max_stars_repo_name: varnerlab/SEML

max_issues_repo_name: varnerlab/SEML

max_forks_repo_name: varnerlab/SEML


text:
using Neighborhood: metricreturntype

const metric = Euclidean()

# Retuns same value as metric but compiler cannot infer type based on arguments
function typeunstable_metric(x, y)
    isnan(x[1]) && return "foo"  # Confuse the compiler
    return metric(x, y)
end

# Wraps a function and records if it was called
struct CalledChecker{F}
    f::F
    called::Ref{Bool}

    CalledChecker(f) = new{typeof(f)}(f, Ref(false))
end

function (c::CalledChecker)(args...)
    c.called[] = true
    return c.f(args...)
end

@testset "metricreturntype" begin
    for T in [Float64, Float32, Int]
        x = zeros(T, 3)

        metric_checked = CalledChecker(metric)
        @test metricreturntype(metric_checked, x) === typeof(metric(x, x))
        @test !metric_checked.called[]  # Function should never have actually been called

        # Return value should itself be inferrable
        @inferred metricreturntype(metric, x)

        # No inference possible, but should still get the correct result by calling the function
        @test metricreturntype(typeunstable_metric, x) === typeof(typeunstable_metric(x, x))
    end
end


text:
# Device type and auxiliary functions

#using .API: VEDAdevice, VEDAdevice_attribute

export
    VEDevice, current_device, has_device, deviceid, physical_deviceid, aveo_deviceid,
    name, totalmem, can_access_peer

"""
    VEDevice(ordinal::Integer)

Get a handle to a compute device.
"""
struct VEDevice
    handle::VEDAdevice

    function VEDevice(ordinal::Integer)
        device_ref = Ref{VEDAdevice}()
        vedaDeviceGet(device_ref, ordinal)
        new(device_ref[])
    end

    global function current_device()
        device_ref = Ref{VEDAdevice}()
        res = vedaCtxGetDevice(device_ref)
        res == VEDA_ERROR_INVALID_CONTEXT && throw(UndefRefError())
        res != VEDA_SUCCESS && throw_api_error(res)
        return _VEDevice(device_ref[])
    end

    # for outer constructors
    global _VEDevice(handle::VEDAdevice) = new(handle)

end

"""
    deviceid(dev::VEDevice)

Returns the ordinal ID of a VEDevice.
"""
function deviceid(dev::VEDevice)
    return dev.handle
end

"""
    physical_deviceid(dev::VEDevice)

Returns the physical device ID of a VEDevice.
"""
function physical_deviceid(dev::VEDevice)
    phys_id = Ref{Int32}()
    vedaDeviceGetPhysicalId(phys_id, dev)
    return phys_id[]
end

"""
    aveo_deviceid(dev::VEDevice)

Returns the AVEO device ID of a VEDevice.
"""
function aveo_deviceid(dev::VEDevice)
    aveo_id = Ref{Int32}()
    vedaDeviceGetAVEOId(aveo_id, dev)
    return aveo_id[]
end

"""
    current_device()

Returns the current device.

!!! warning

    This is a low-level API, returning the current device as known to the CUDA driver.
    For most users, it is recommended to use the [`device`](@ref) method instead.
"""
current_device()

"""
    has_device()

Returns whether there is an active device.
"""
function has_device()
    device_ref = Ref{VEDAdevice}()
    res = vedaCtxGetDevice(device_ref)
    if res == VEDA_SUCCESS
        return true
    elseif res == VEDA_ERROR_INVALID_CONTEXT
        return false
    else
        throw_api_error(res)
    end
end

Base.convert(::Type{VEDAdevice}, dev::VEDevice) = dev.handle

function Base.show(io::IO, ::MIME"text/plain", dev::VEDevice)
    print(io, "VEDevice($(dev.handle)): ")
    print(io, "physical_id=$(physical_deviceid(dev)) ")
    print(io, "aveo_id=$(aveo_deviceid(dev)) ")
    print(io, "$(name(dev))")
end

"""
    name(dev::VEDevice)

Returns an identifier string for the device.
"""
function name(dev::VEDevice)
    buflen = 256
    buf = Vector{Cchar}(undef, buflen)
    vedaDeviceGetName(pointer(buf), buflen, dev)
    buf[end] = 0
    return unsafe_string(pointer(buf))
end

"""
    totalmem(dev::VEDevice)

Returns the total amount of memory (in bytes) on the device.
"""
function totalmem(dev::VEDevice)
    mem_ref = Ref{Csize_t}()
    vedaDeviceTotalMem_v2(mem_ref, dev)
    return mem_ref[]
end

# TODO: in principle this is true, but quick transfers are not yet
#       in AVEO.
function can_access_peer(dev::VEDevice, peer::VEDevice)
    return false
end

## device iteration

export devices, ndevices

struct DeviceIterator end

"""
    devices()

Get an iterator for the compute devices.
"""
devices() = DeviceIterator()

Base.eltype(::DeviceIterator) = VEDevice

function Base.iterate(iter::DeviceIterator, i=1)
    i >= length(iter) + 1 ? nothing : (VEDevice(i-1), i+1)
end

Base.length(::DeviceIterator) = ndevices()

Base.IteratorSize(::DeviceIterator) = Base.HasLength()

function Base.show(io::IO, ::MIME"text/plain", iter::DeviceIterator)
    print(io, "VEDA.DeviceIterator() for $(length(iter)) devices:\n")
    if !isempty(iter)
        for dev in iter
            print(io, "$(deviceid(dev)). ")
            print(io, "physical_id=$(physical_deviceid(dev)) ")
            print(io, "aveo_id=$(aveo_deviceid(dev)) ")
            print(io, "$(name(dev))\n")
        end
    end
end

function ndevices()
    count_ref = Ref{Cint}()
    vedaDeviceGetCount(count_ref)
    return count_ref[]
end

## attributes

export attribute

"""
    attribute(dev::VEDevice, code)

Returns information about the device.
"""
function attribute(dev::VEDevice, code::VEDAdevice_attribute)
    value_ref = Ref{Cint}()
    vedaDeviceGetAttribute(value_ref, code, dev)
    return value_ref[]
end

@enum_without_prefix VEDAdevice_attribute VEDA_DEVICE_ATTRIBUTE_


text:
export WeightedMeasure

"""
    struct WeightedMeasure{R,M} <: AbstractMeasure
        logweight :: R
        base :: M
    end
"""

abstract type AbstractWeightedMeasure <: AbstractMeasure end

logweight(μ::AbstractWeightedMeasure) = μ.logweight
basemeasure(μ::AbstractWeightedMeasure) = μ.base

TV.as(μ::AbstractWeightedMeasure) = TV.as(μ.base)

function logdensity(sm::AbstractWeightedMeasure, x)
    logdensity(sm.base, x) + sm.logweight
end

###############################################################################

struct WeightedMeasure{R,M} <: AbstractWeightedMeasure
    logweight :: R
    base :: M
end

function Base.show(io::IO, μ::WeightedMeasure)
    io = IOContext(io, :compact => true)
    print(io, exp(μ.logweight), " * ", μ.base)
end

function Base.show_unquoted(io::IO, μ::WeightedMeasure, indent::Int, prec::Int)
    io = IOContext(io, :compact => true)
    if Base.operator_precedence(:*) ≤ prec
        print(io, "(")
        show(io, μ)
        print(io, ")")
    else
        show(io, μ)
    end
    return nothing
end

function Base.:*(k::T, m::AbstractMeasure) where {T <: Number}
    logk = log(k)
    return WeightedMeasure{typeof(logk), typeof(m)}(logk,m)
end

Base.:*(m::AbstractMeasure, k::Real) = k * m

≪(::M, ::WeightedMeasure{R,M}) where {R,M} = true
≪(::WeightedMeasure{R,M}, ::M) where {R,M} = true

sampletype(μ:: WeightedMeasure) = sampletype(μ.base)

###############################################################################
struct ParamWeightedMeasure{F,N,T,R,B} <: AbstractWeightedMeasure
    f::F
    par::NamedTuple{N,T}
    logweight::R
    base::B

    function ParamWeightedMeasure(f::F, par::NamedTuple{N,T}, base::B) where {F,N,T,B}
        ℓ = f(par)
        new{F,N,T,typeof(ℓ),B}(f,par,ℓ,base)
    end
end

function Base.show(io::IO, d::ParamWeightedMeasure)
    io = IOContext(io, :compact => true)
    print(io, "ParamWeighted(",d.f, ", ", d.par,", ", d.base, ")")
end


text:
function LinearAlgebra.normalize!(sim::Simulation{1},ψ,B)
    n = length(sim.domain_indices)
    lattices = map(d->d.lattice,sim.domains)
    dx = map(i->lattices[i].dx,sim.domain_indices)
    for i ∈ 1:size(ψ,2)
        𝒩² = zero(eltype(ψ))
        for j ∈ 1:size(ψ,1)
            𝒩² += (ψ[j,i].^2)*B[j,j]*dx[mod1(j,n)]
        end
        ψ[:,i] /= sqrt(𝒩²)*exp(complex(0,angle(ψ[end÷2-1,i])))
    end
    return nothing
end

"""
    eig_kl(sim, ω; ky=0, kz=0, verbose=false, kwargs...]) -> ω, ψ

Linear eigenvalue solve for `ω`.

Keyword `verbose` defaults to `false`.
See docs for ArnoldiMethodWrapper for details of `kwargs`.
"""
function eig_kl(sim::Simulation{1,Symmetric,T},
            ω::Number;
            ky::Number=0,
            kz::Number=0,
            verbose::Bool=false,
            kwargs...
            ) where T

    A,B = maxwell_lep(sim;ky=ky,kz=kz)(ω)
    decomp, history = partialschur(A, B, ω^2; diag_inv_B=true, kwargs...)
    history.converged || @warn "$(history.nev - history.nconverged) did not converge"
    verbose ? println(history) : nothing

    λ::Array{T,1}, v::Array{T,2} = partialeigen(decomp, ω^2)
    normalize!(sim,v,B) # Normalize according to (ψ₁,ψ₂)=δ₁₂
    return sqrt.(λ), v
end

"""
    eig_knl(sim, k, ka=0, kb=0; method=contour_beyn, nk=3, display=false, quad_n=100, kwargs...) -> k,ψ
"""
function eig_knl(
            sim::Simulation{1,Symmetric,T},
            ω::Number;
            ky::Number=0,
            kz::Number=0,
            verbose::Bool = false,
            check_consistency::Bool = true,
            Schur_fact::Bool = false,
            align_sparsity_patterns::Bool = false,
            nlmethod::Function = iar,
            logger::Integer = Int(verbose),
            kwargs...
            ) where T

    spmf_kwargs = (
        :check_consistency=>check_consistency,
        :Schur_fact=>Schur_fact,
        :align_sparsity_patterns=>align_sparsity_patterns)
    nep = maxwell_nep(sim;
        ky=ky,
        kz=kz,
        :check_consistency => check_consistency,
        :Schur_fact => Schur_fact,
        :align_sparsity_patterns => align_sparsity_patterns)

    λ::Array{T,1}, v::Array{T,2} = nlmethod(nep; σ=ω, kwargs...,logger=logger)
    return λ, v
end

"""
    eig_cf(sim, k, [ka=0, kb=0; η=0, verbose, lupack, kwargs...]) -> η, u

Linear CF eigenvalues closest to `η`

Keyword `verbose` defaults to `false`.
Keyword `lupack` defaults to `:auto` and contrls which package is used in the eigensolver.
See docs for ArnoldiMethodWrapper for details.
"""
function eig_cf(sim::Simulation{1,Symmetric,T},
            ω::Number;
            ky::Number=0,
            kz::Number=0,
            η::Number=0,
            verbose::Bool=false,
            kwargs...
            ) where T

    A,B = maxwell_lep(sim)(ω)
    A = A-B*ω^2
    F = -spdiagm(0=>vcat(sim.α[1].*sim.F,sim.α[1].*sim.F,sim.α[1].*sim.F))
    decomp, history = partialschur(A, F*ω^2, η; diag_inv_B=false, kwargs...)
    history.converged || @warn "$(history.nev-history.nconverged) eigenvectors failed to converge"
    verbose ? println(history) : nothing

    λ::Array{T,1}, v::Array{T,2} = partialeigen(decomp,η)
    normalize!(sim,v,B)
    return λ, v
end


text:
count(vs::VarStats) = vs.n_obs
sum(vs::VarStats) = vs.sum
min(vs::VarStats) = vs.min
max(vs::VarStats) = vs.max
mean(vs::VarStats) = vs.sum / vs.n_obs
var(vs::VarStats) = vs.var
std(vs::VarStats) = sqrt(vs.var)

function updatestats!(vs::VarStats, value::Real)
    vs.n_obs += 1
    vs.sum += value
    if value < vs.min
        vs.min = value
    end
    if value > vs.max
        vs.max = value
    end
    if vs.n_obs > 2
        n = vs.n_obs
        vs.var = (n - 2) / (n - 1) * vs.n_obs + 1 / n * (value - mean(vs))^2
    end
end

function updatestats!{T<:Real}(varstats::Vector{VarStats}, values::Vector{T})
    for i=1:length(varstats)
        updatestats!(varstats[i], values[i])
    end
end

function aggregate!{N<:Number}(agg::StatsAggregator,
                               key::Any, values::Vector{N})
    if !haskey(agg.groups, key)
        # init new varstats group to the length of first incoming values
        agg.groups[key] = [VarStats() for i=1:length(values)]
    end
    updatestats!(agg.groups[key], values)
end

haskey(agg::StatsAggregator, key::Any) = haskey(agg.groups, key)
keys(agg::StatsAggregator) = keys(agg.groups)
values(agg::StatsAggregator) = values(agg.groups)
getindex(agg::StatsAggregator, k::Any) = getindex(agg.groups, k)


text:
using LinearAlgebra: I
using TailRiskScenGen
using Distributions

function valid_scenario_set(scenarios::Matrix{Float64}, probs::Vector{Float64})
    @test size(scenarios, 2) == length(probs)
    @test sum(probs) ≈ 1.0
    for p in probs
        @test (p >= 0)
    end
end

d = 10
num_scen = 1000
μ = rand(Normal(), d)
A = rand(d,d)./d
Σ = A'A
K = FiniteCone(-Array{Float64}(I, d, d))
dist = MvNormal(μ, Σ)
Ω = EllipticalRiskRegion(dist, K, 0.95)
scenarios = rand(dist, num_scen)

new_scen, new_prob = aggregate_scenarios(scenarios, Ω)
valid_scenario_set(new_scen, new_prob)

agg_sample_scen, agg_sample_prob = aggregation_sampling(dist, Ω, num_scen)
valid_scenario_set(agg_sample_scen, agg_sample_prob)

cluster_sample_scen, cluster_sample_prob = nonrisk_clustering(dist, Ω, 100,10)
valid_scenario_set(cluster_sample_scen, cluster_sample_prob)

cluster_scen, cluster_probs = nonrisk_clustering(scenarios, Ω, 10)
valid_scenario_set(cluster_scen, cluster_probs)

# Monotonic Risk Region
sample = rand(dist, 10000)
α = 0.99
surv = SurvivorApproximator(sample, α)
Ω = MonotonicRiskRegion(surv, 0.99)
new_scen, new_prob = aggregate_scenarios(scenarios, Ω)
valid_scenario_set(new_scen, new_prob)

agg_sample_scen, agg_sample_prob = aggregation_sampling(dist, Ω, num_scen)
valid_scenario_set(agg_sample_scen, agg_sample_prob)

cluster_sample_scen, cluster_sample_prob = nonrisk_clustering(dist, Ω, 100,10)
valid_scenario_set(cluster_sample_scen, cluster_sample_prob)

cluster_scen, cluster_probs = nonrisk_clustering(scenarios, Ω, 10)
valid_scenario_set(cluster_scen, cluster_probs)


text:
module PyThermo

using PyCall
using Printf
using Unitful

import Base: convert, ==, isequal, hash, getindex, setindex!, haskey, keys, show

export Thermo, ShockTube
export Species, Mixture
export isentropic_exponent, temperature, pressure, density, molar_density, R_specific, soundspeed

include("init.jl")
const Thermo = PY_THERMO

_unit(x, u) = Float64(ustrip(uconvert(u, x)))

function _SI_TP(kwargs)
    d = Dict{Any, Any}(kwargs)
    T = get!(d, :T, 298.15)
    P = get!(d, :P, 101325.0)
    d[:T] = (T isa Unitful.Temperature) ? _unit(T, u"K") : Float64(T)
    d[:P] = (P isa Unitful.Pressure)    ? _unit(P, u"Pa") : Float64(P)
    return (d...,)
end

abstract type Chemical end

"""
    Species(ID::String; T=298.15, P=101325)

Creates a `Species` object which contains basic information such as
molecular weight and the structure of the species, as well as thermodynamic
and transport properties as a function of temperature and pressure.
If `T` and `P` are given as non-`Unitful` numbers, they must have units of `K` and `Pa`.

Parameters
----------
ID : One of the following [-]:

* Name, in IUPAC form or common form or a synonym registered in PubChem \n
* InChI name, prefixed by "InChI=1S/" or "InChI=1/" \n
* InChI key, prefixed by "InChIKey=" \n
* PubChem CID, prefixed by "PubChem=" \n
* SMILES (prefix with "SMILES=" to ensure smiles parsing) \n
* CAS number

T : temperature of the chemical (default 298.15 K) \n
P : pressure of the chemical (default 101325 Pa)

Examples
--------
```jldoctest
julia> He = Species("He")
Species(He, 298.1 K, 1.013e+05 Pa)

julia> density(He)
0.16352343013746262 kg m^-3

julia> using Unitful

julia> He.T = 30u"K"
30 K

julia> density(He)
1.6235030074934973 kg m^-3
```

A wide variety of unexported properties can be accessed from the underlying Python object:
```
julia> SF6 = Species("SF6", P=30u"psi")
Species(SF6, 298.1 K, 2.068e+05 Pa)

julia> SF6.MW
146.055419

help?> SF6.<tab>
A                                 SGs                                __init_subclass__
API                               STEL                               __le__
Am                                STEL_source                        __lt__
Bond                              STEL_sources                       __module__
Bvirial                           S_dep_Tb_P_ref_g                   __ne__
CAS                               S_dep_Tb_Pb_g                      __new__
Capillary                         S_dep_Tb_Pb_l                      __reduce__
Carcinogen                        S_dep_ref_g                        __reduce_ex__
Carcinogen_source                 S_int_Tb_to_T_ref_g                __repr__
Carcinogen_sources                S_int_l_Tm_to_Tb                   __setattr__
[...]
```
"""
struct Species <: Chemical
    o::PyObject
end
Species(chemname::String; kwargs...) = Species(PY_CHEM.Chemical(chemname; _SI_TP(kwargs)...))

Base.show(io::IO, species::Species) = @printf(io, "Species(%s, %0.1f K, %0.3e Pa)", species.ID, species.T, species.P)

"""
    Mixture(chemnames::Vector{String}; kwargs...)
    Mixture(chemnames::Vector{Pair{String, Float64}}, kwargs...)

Creates a `Mixture` object which contains basic information such as
molecular weight and the structure of the species, as well as thermodynamic
and transport properties as a function of temperature and pressure.

The components of the mixture are specified by the names of
the chemicals; the composition can be specified by providing any one of the
following parameters as a keyword argument:

* Mass fractions `ws`
* Mole fractions `zs`
* Liquid volume fractions (based on pure component densities) `Vfls`
* Gas volume fractions (based on pure component densities) `Vfgs`

The composition can also be specified by providing a vector of `"ID" => molefrac` pairs.

Examples
--------
```jldoctest
julia> air = Mixture(["N2" => 0.78, "O2" => 0.21, "Ar" => 0.01])
Mixture({N2: 0.78, O2: 0.21, Ar: 0.01}, 298.1 K, 1.013e+05 Pa)

julia> soundspeed(air)
346.1466532754559 m s^-1
```
"""
struct Mixture <: Chemical
    o::PyObject
end
Mixture(chemnames::Vector{String}; kwargs...) = Mixture(PY_CHEM.Mixture(chemnames; _SI_TP(kwargs)...))

function Mixture(chems::Vector{Pair{String, Float64}}; kwargs...)
    Mixture(first.(chems); kwargs..., zs = last.(chems))
end

function Mixture(chems::Vector{Pair{T1, T2}}; kwargs...) where {T1, T2}
    Mixture(string.(first.(chems)) .=> Float64.(last.(chems)); kwargs...)
end

# Mixture(args...; kwargs...) = Mixture([args...], kwargs...)

function composition_string(mix)
    species_str = try
        mix.IDs
    catch
        mix.names
    end
    s = "{"
    for (species, χ) in zip(species_str, mix.zs)
        s *= @sprintf("%s: %0.3g, ", species, χ)
    end
    s[1:end-2] * "}"
end
Base.show(io::IO, mix::Mixture) = @printf(io, "Mixture(%s, %0.1f K, %0.3e Pa)", composition_string(mix), mix.T, mix.P)

PyObject(c::Chemical) = getfield(c, :o)
convert(::Type{T}, o::PyCall.PyObject) where {T <: Chemical} = T(o)
==(c1::Chemical, c2::Chemical) = PyObject(c1) == PyObject(c2)
hash(c::Chemical) = hash(PyObject(c))
Base.copy(c::T) where {T <: Chemical} = T(PY_COPY.copy(PyObject(c)))

Base.Docs.doc(c::Chemical) = Base.Docs.doc(PyObject(c))
Base.Docs.Binding(c::Chemical, s::Symbol) = Base.Docs.Binding(PyObject(c), s)

Base.getproperty(c::Chemical, s::Symbol) = getproperty(PyObject(c), s)
Base.setproperty!(c::Chemical, s::Symbol, x) = setproperty!(PyObject(c), s, x)
Base.hasproperty(c::Chemical, s::Symbol) = hasproperty(PyObject(c), s)
Base.propertynames(c::Chemical) = propertynames(PyObject(c))
haskey(c::Chemical, x) = haskey(PyObject(c), x)

# Strip units for unitful `setproperty`
Base.setproperty!(c::Chemical, s::Symbol, T::Unitful.Temperature) = setproperty!(PyObject(c), s, _unit(T, u"K"))
Base.setproperty!(c::Chemical, s::Symbol, T::Unitful.Pressure) = setproperty!(PyObject(c), s, _unit(T, u"Pa"))

# Thermodynamic property accessors
temperature(c::Chemical)  = c.T * u"K"
pressure(c::Chemical)  = c.P * u"Pa"
density(c::Chemical)  = c.rho * u"kg/m^3"
molar_density(c::Chemical) = c.rhom * u"mol/m^3"
isentropic_exponent(c::Chemical)  = c.isentropic_exponent
R_specific(c::Chemical)  = c.R_specific * u"J/(kg*K)"
soundspeed(c::Chemical) = sqrt(isentropic_exponent(c) * R_specific(c) * temperature(c)) |> u"m/s"

include("ShockTube.jl")
end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

@with_kw mutable struct ConstantStepParameters{T <: AbstractFloat}
    γ::T = 0.01
end

"""
    ConstantStep

Functor object for using a constant step size in a quasigradient algorithm. Create by supplying a [`Constant`](@ref) object through `step` to `QuasiGradient.Optimizer` or by setting the [`StepSize`](@ref) attribute.

...
# Parameters
- `γ::AbstractFloat = 0.01`: Step length
...
"""
struct ConstantStep{T <: AbstractFloat} <: AbstractStep
    parameters::ConstantStepParameters{T}

    function ConstantStep(::Type{T}; kw...) where T <: AbstractFloat
        return new{T}(ConstantStepParameters{T}(; kw...))
    end
end
step(step::ConstantStep, ::Integer, ::Float64, ::AbstractVector, ::AbstractVector) = step.parameters.γ

# API
# ------------------------------------------------------------
"""
    Constant

Factory object for [`ConstantStep`](@ref). Pass to `step` in `Quasigradient.Optimizer` or set the [`StepSize`](@ref) attribute. See ?ConstantStep for parameter descriptions.

"""
struct Constant <: AbstractStepSize
    parameters::ConstantStepParameters{Float64}
end
Constant(γ::AbstractFloat) = Constant(ConstantStepParameters(; γ = Float64(γ)))
Constant(; kw...) = Constant(ConstantStepParameters(; kw...))

function (stepsize::Constant)(::Type{T}) where T <: AbstractFloat
    return ConstantStep(T; type2dict(stepsize.parameters)...)
end

function str(::Constant)
    return "constant step size"
end


text:
export VBasedPolicy

"""
    VBasedPolicy(learner::V, f::F) -> VBasedPolicy{V,F}

Use a function `f` to generate action and update the V-`learner` in the meantime.
"""
Base.@kwdef struct VBasedPolicy{V,F} <: AbstractPolicy
    learner::V
    f::F
end

(π::VBasedPolicy)(obs) = π.f(obs)

get_prob(π::VBasedPolicy, s, a) = get_prob(π.f, s, a)

extract_transitions(buffer, π::VBasedPolicy) = extract_transitions(buffer, π.learner)

function extract_transitions(
    buffer::EpisodeTurnBuffer,
    ::MonteCarloLearner{T,A},
) where {T,A<:AbstractVApproximator}
    if isfull(buffer)
        @views (states = state(buffer)[1:end-1], rewards = reward(buffer)[2:end])
    else
        nothing
    end
end