text:
#---------------------------------------------------------------------#
#This function computes the high-order grid & elements in 3D.
#Written by F.X. Giraldo on 5/2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#

using Plots
include("warp_grid.jl")

function create_grid(Np,Npoin,Ne,Nboun,Nex,Ney,Nez,ξ,plot_grid,warp_grid,DFloat)

    #Initialize Global Arrays
    coord=zeros(DFloat,3,Npoin)
    intma=zeros(Int64,Np,Np,Np,Ne)
    iboun=zeros(Int64,Nboun,1)

    #Set some constants
    xmin=-1
    xmax=+1
    ymin=-1
    ymax=+1
    zmin=-1
    zmax=+1
    dx=(xmax-xmin)/Nex
    dy=(ymax-ymin)/Ney
    dz=(zmax-zmin)/Nez
    N=Np-1
    Nx=Nex*N + 1
    Ny=Ney*N + 1
    Nz=Nez*N + 1

    #Initialize Local Arrays
    node=zeros(Int64,Nx,Ny,Nz)

    #GENERATE COORD
    I=0
    kk=0
    for m=1:Nez
        z0=zmin + (m-1)*dz
        if (m == 1)
            n1=1
        else
            n1=2
        end
        for n=n1:Np
            z=( ξ[n]+1 )*dz/2 + z0
            kk=kk+1
            jj=0
            for k=1:Ney
                y0=ymin + (k-1)*dy
                if (k == 1)
                    l1=1
                else
                    l1=2
                end
                for l=l1:Np
                    y=( ξ[l]+1 )*dy/2 + y0
                    jj=jj+1
                    ii=0
                    for i=1:Nex
                        x0=xmin + (i-1)*dx
                        if (i == 1)
                            j1=1
                        else
                            j1=2
                        end
                        for j=j1:Np
                            ii=ii + 1
                            I=I + 1
                            x=( ξ[j]+1 )*dx/2 + x0
                            coord[1,I]=x
                            coord[2,I]=y
                            coord[3,I]=z
                            node[ii,jj,kk]=I
                        end #j
                    end #i
                end #l
            end #k
        end #m
    end #n

    #GENERATE INTMA
    e=0
    for m=1:Nez
        for k=1:Ney
            for i=1:Nex
                e=e+1
                for n=1:Np
                    kk=N*(m-1) + n
                    for l=1:Np
                        jj=N*(k-1) + l
                        for j=1:Np
                            ii=N*(i-1) + j
                            I=node[ii,jj,kk]
                            intma[j,l,n,e]=I
                        end #j
                    end #l
                end #n
            end #i
        end #k
    end #m

    #GENERATE Boundary Conditions
    B=0

    #Front Face [y=-1]
    j=1
    for k=1:Nz
        for i=1:Nx
            I=node[i,j,k]
            B=B + 1
            iboun[B]=I
        end
    end

    #Back Face [y=+1]
    j=Ny
    for k=1:Nz
        for i=1:Nx
            I=node[i,j,k]
            B=B + 1
            iboun[B]=I
        end
    end

    #Left Face [x=-1]
    i=1
    for k=1:Nz
        for j=2:Ny-1
            I=node[i,j,k]
            B=B + 1
            iboun[B]=I
        end
    end

    #Right Face [x=+1]
    i=Nx
    for k=1:Nz
        for j=2:Ny-1
            I=node[i,j,k]
            B=B + 1
            iboun[B]=I
        end
    end

    #Bottom Face [z=-1]
    k=1
    for j=2:Ny-1
        for i=2:Nx-1
            I=node[i,j,k]
            B=B + 1
            iboun[B]=I
        end
    end

    #Top Face [z=-1]
    k=Nz
    for j=2:Ny-1
        for i=2:Nx-1
            I=node[i,j,k]
            B=B + 1
            iboun[B]=I
        end
    end

    if (B != Nboun)
        disp("Error in Create_Grid")
        return
    end

    #Warp Grid
    if (warp_grid)
        warp_grid!(coord,Npoin,DFloat)
    end

    #Plot Grid
    title_text=["Grid : Ne =  $Ne,  N = $N"]
    if (plot_grid)
        x=zeros(DFloat,5,1)
        y=zeros(DFloat,5,1)
        z=zeros(DFloat,5,1)
        x[1]=-1; y[1]=-1; z[1]=-1
        x[2]=+1; y[2]=-1; z[2]=-1
        x[3]=+1; y[3]=+1; z[3]=-1
        x[4]=-1; y[4]=+1; z[4]=-1
        x[5]=-1; y[5]=-1; z[5]=-1
        p1=plot(x,y,z,linecolor="black",legend=false)

        for e=1:Ne

            #Build x-y plane grid for each z value
            for k=1:Np
                for j=1:N
                    for i=1:N
                        i1=intma[i,j,k,e]
                        i2=intma[i+1,j,k,e]
                        i3=intma[i+1,j+1,k,e]
                        i4=intma[i,j+1,k,e]
                        x[1]=coord[1,i1]; y[1]=coord[2,i1]; z[1]=coord[3,i1]
                        x[2]=coord[1,i2]; y[2]=coord[2,i2]; z[2]=coord[3,i2]
                        x[3]=coord[1,i3]; y[3]=coord[2,i3]; z[3]=coord[3,i3]
                        x[4]=coord[1,i4]; y[4]=coord[2,i4]; z[4]=coord[3,i4]
                        x[5]=coord[1,i1]; y[5]=coord[2,i1]; z[5]=coord[3,i1]
                        p1=plot!(p1,x,y,z,linecolor="red",legend=false)
                    end
                end
            end

            #Build x-z plane grid for each y value
            for j=1:Np
                for k=1:N
                    for i=1:N
                        i1=intma[i,j,k,e]
                        i2=intma[i+1,j,k,e]
                        i3=intma[i+1,j,k+1,e]
                        i4=intma[i,j,k+1,e]
                        x[1]=coord[1,i1]; y[1]=coord[2,i1]; z[1]=coord[3,i1]
                        x[2]=coord[1,i2]; y[2]=coord[2,i2]; z[2]=coord[3,i2]
                        x[3]=coord[1,i3]; y[3]=coord[2,i3]; z[3]=coord[3,i3]
                        x[4]=coord[1,i4]; y[4]=coord[2,i4]; z[4]=coord[3,i4]
                        x[5]=coord[1,i1]; y[5]=coord[2,i1]; z[5]=coord[3,i1]
                        p1=plot!(p1,x,y,z,linecolor="red",legend=false)
                    end
                end
            end

            #Build y-z plane grid for each x value
            for i=1:Np
                for k=1:N
                    for j=1:N
                        i1=intma[i,j,k,e]
                        i2=intma[i,j+1,k,e]
                        i3=intma[i,j+1,k+1,e]
                        i4=intma[i,j,k+1,e]
                        x[1]=coord[1,i1]; y[1]=coord[2,i1]; z[1]=coord[3,i1]
                        x[2]=coord[1,i2]; y[2]=coord[2,i2]; z[2]=coord[3,i2]
                        x[3]=coord[1,i3]; y[3]=coord[2,i3]; z[3]=coord[3,i3]
                        x[4]=coord[1,i4]; y[4]=coord[2,i4]; z[4]=coord[3,i4]
                        x[5]=coord[1,i1]; y[5]=coord[2,i1]; z[5]=coord[3,i1]
                        p1=plot!(p1,x,y,z,linecolor="red",legend=false)
                   end
                end
            end

            #---------------------Build Element Boundary----------------------#
            #Build x-y plane grid for k=1 value
            i1=intma[1,1,1,e]
            i2=intma[Np,1,1,e]
            i3=intma[Np,Np,1,e]
            i4=intma[1,Np,1,e]
            x[1]=coord[1,i1]; y[1]=coord[2,i1]; z[1]=coord[3,i1]
            x[2]=coord[1,i2]; y[2]=coord[2,i2]; z[2]=coord[3,i2]
            x[3]=coord[1,i3]; y[3]=coord[2,i3]; z[3]=coord[3,i3]
            x[4]=coord[1,i4]; y[4]=coord[2,i4]; z[4]=coord[3,i4]
            x[5]=coord[1,i1]; y[5]=coord[2,i1]; z[5]=coord[3,i1]
            p1=plot!(p1,x,y,z,linecolor="blue",legend=false)

            #Build x-y plane grid for k=Np value
            i1=intma[1,1,Np,e]
            i2=intma[Np,1,Np,e]
            i3=intma[Np,Np,Np,e]
            i4=intma[1,Np,Np,e]
            x[1]=coord[1,i1]; y[1]=coord[2,i1]; z[1]=coord[3,i1]
            x[2]=coord[1,i2]; y[2]=coord[2,i2]; z[2]=coord[3,i2]
            x[3]=coord[1,i3]; y[3]=coord[2,i3]; z[3]=coord[3,i3]
            x[4]=coord[1,i4]; y[4]=coord[2,i4]; z[4]=coord[3,i4]
            x[5]=coord[1,i1]; y[5]=coord[2,i1]; z[5]=coord[3,i1]
            p1=plot!(p1,x,y,z,linecolor="blue",legend=false)

            #Build x-z plane grid for j=1
            i1=intma[1,1,1,e]
            i2=intma[Np,1,1,e]
            i3=intma[Np,1,Np,e]
            i4=intma[1,1,Np,e]
            x[1]=coord[1,i1]; y[1]=coord[2,i1]; z[1]=coord[3,i1]
            x[2]=coord[1,i2]; y[2]=coord[2,i2]; z[2]=coord[3,i2]
            x[3]=coord[1,i3]; y[3]=coord[2,i3]; z[3]=coord[3,i3]
            x[4]=coord[1,i4]; y[4]=coord[2,i4]; z[4]=coord[3,i4]
            x[5]=coord[1,i1]; y[5]=coord[2,i1]; z[5]=coord[3,i1]
            p1=plot!(p1,x,y,z,linecolor="blue",legend=false)

            #Build x-z plane grid for j=Np
            i1=intma[1,Np,1,e]
            i2=intma[Np,Np,1,e]
            i3=intma[Np,Np,Np,e]
            i4=intma[1,Np,Np,e]
            x[1]=coord[1,i1]; y[1]=coord[2,i1]; z[1]=coord[3,i1]
            x[2]=coord[1,i2]; y[2]=coord[2,i2]; z[2]=coord[3,i2]
            x[3]=coord[1,i3]; y[3]=coord[2,i3]; z[3]=coord[3,i3]
            x[4]=coord[1,i4]; y[4]=coord[2,i4]; z[4]=coord[3,i4]
            x[5]=coord[1,i1]; y[5]=coord[2,i1]; z[5]=coord[3,i1]
            p1=plot!(p1,x,y,z,linecolor="blue",legend=false)

            #Build y-z plane grid for each i=1 value
            i1=intma[1,1,1,e]
            i2=intma[1,Np,1,e]
            i3=intma[1,Np,Np,e]
            i4=intma[1,1,Np,e]
            x[1]=coord[1,i1]; y[1]=coord[2,i1]; z[1]=coord[3,i1]
            x[2]=coord[1,i2]; y[2]=coord[2,i2]; z[2]=coord[3,i2]
            x[3]=coord[1,i3]; y[3]=coord[2,i3]; z[3]=coord[3,i3]
            x[4]=coord[1,i4]; y[4]=coord[2,i4]; z[4]=coord[3,i4]
            x[5]=coord[1,i1]; y[5]=coord[2,i1]; z[5]=coord[3,i1]
            p1=plot!(p1,x,y,z,linecolor="blue",legend=false)

            #Build y-z plane grid for each i=Np value
            i1=intma[Np,1,1,e]
            i2=intma[Np,Np,1,e]
            i3=intma[Np,Np,Np,e]
            i4=intma[Np,1,Np,e]
            x[1]=coord[1,i1]; y[1]=coord[2,i1]; z[1]=coord[3,i1]
            x[2]=coord[1,i2]; y[2]=coord[2,i2]; z[2]=coord[3,i2]
            x[3]=coord[1,i3]; y[3]=coord[2,i3]; z[3]=coord[3,i3]
            x[4]=coord[1,i4]; y[4]=coord[2,i4]; z[4]=coord[3,i4]
            x[5]=coord[1,i1]; y[5]=coord[2,i1]; z[5]=coord[3,i1]
            p1=plot!(p1,x,y,z,linecolor="blue",legend=false)
        end #e
        display(plot(p1))
    end #plot_grid

    return(coord,intma,iboun)
end


text:
# core parts of add, sub, mul, div

@inline function addovf(x, y)
    ovf = false
    numer, ovfl = mul_with_overflow(x.num, x.den) # here, numer is a temp
    ovf |= ovfl
    denom, ovfl = mul_with_overflow(x.den, x.num) # here, denom is a temp
    ovf |= ovfl
    numer, ovfl = add_with_overflow(numer, denom) # numerator of sum
    ovf |= ovfl
    denom, ovfl = mul_with_overflow(x.den, y.den) # denominator of sum
    ovf |= ovfl
    return numer, denom, ovf
end

@inline function subovf(x, y)
    ovf = false
    numer, ovfl = mul_with_overflow(x.num, x.den) # here, numer is a temp
    ovf |= ovfl
    denom, ovfl = mul_with_overflow(x.den, y.num) # here, denom is a temp
    ovf |= ovfl
    numer, ovfl = sub_with_overflow(numer, denom) # numerator of difference
    ovf |= ovfl
    denom, ovfl = mul_with_overflow(x.den, y.den) # denominator of difference
    ovf |= ovfl
    return numer, denom, ovf
end

@inline function mulovf(x, y)
    ovf = false
    numer, ovfl = mul_with_overflow(x.num, y.num)
    ovf |= ovfl
    denom, ovfl = mul_with_overflow(x.den, y.den)
    ovf |= ovfl
    return numer, denom, ovf
end

@inline function divovf(x, y)
    ovf = false
    numer, ovfl = mul_with_overflow(x.num, y.den)
    ovf |= ovfl
    denom, ovfl = mul_with_overflow(x.den, y.num)
    ovf |= ovfl
    return numer, denom, ovf
end

ext: jl


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -12.75118444054646 <= q <= 9.039664285972906)

                     add_NL_constraint(m, :(sigmoid(0.715796138346545 + -0.1516695982843861*sigmoid(-0.317423658088571 + -0.8843064794993589*sigmoid(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + 0.6219431368805965*sigmoid(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + -0.24757066957523488*sigmoid(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))) + -0.0820632819191216*sigmoid(-0.8960501996212993 + -0.7682728024607406*sigmoid(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + -0.3378399703859083*sigmoid(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + 0.4153867033969405*sigmoid(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))) + -0.0571486344530765*sigmoid(0.054021626349418206 + -0.05819989275612114*sigmoid(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + 0.8875139231924525*sigmoid(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + 0.11280184637312685*sigmoid(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3])))) + sigmoid(-0.8464220616292248 + -0.9331432460947604*sigmoid(-0.317423658088571 + -0.8843064794993589*sigmoid(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + 0.6219431368805965*sigmoid(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + -0.24757066957523488*sigmoid(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))) + 0.425513343659472*sigmoid(-0.8960501996212993 + -0.7682728024607406*sigmoid(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + -0.3378399703859083*sigmoid(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + 0.4153867033969405*sigmoid(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))) + 0.1421259709107252*sigmoid(0.054021626349418206 + -0.05819989275612114*sigmoid(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + 0.8875139231924525*sigmoid(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + 0.11280184637312685*sigmoid(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3])))) + sigmoid(-0.5205460270098605 + 0.13812073353255405*sigmoid(-0.317423658088571 + -0.8843064794993589*sigmoid(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + 0.6219431368805965*sigmoid(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + -0.24757066957523488*sigmoid(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))) + 0.7749706249828683*sigmoid(-0.8960501996212993 + -0.7682728024607406*sigmoid(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + -0.3378399703859083*sigmoid(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + 0.4153867033969405*sigmoid(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))) + 0.5007112790542019*sigmoid(0.054021626349418206 + -0.05819989275612114*sigmoid(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + 0.8875139231924525*sigmoid(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + 0.11280184637312685*sigmoid(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
abstract SalvagingAlgo
type NoneSalvagingAlgo <: SalvagingAlgo end

function normalize_pseudo_root!(matrix::Matrix, pseudo::Matrix)
  sz = size(matrix)[1]
  sz == size(pseudo)[1] || error("matrix / pseudo mismatch")

  pseudoCols = size(pseudo)[2]
  # row normalization
  for i = 1:sz
    norm = 0.0
    for j = 1:pseudoCols
      norm += pseudo[i,j] * pseudo[i, j]
    end

    if norm > 0.0
      normAdj = sqrt(matrix[i, i] / norm)
      for j = 1:pseudoCols
        pseudo[i, j] *= normAdj
      end
    end
  end

  return pseudo
end

function rank_reduced_sqrt(matrix::Matrix, maxRank::Int, componentRetainedPercentage::Float64, sa::NoneSalvagingAlgo)
  sz = size(matrix)[1]
  # Symmetric Schur decomp
  jd = schurfact(matrix)
  eigenVals = jd[:values]

  # eigen vals are sorted in decreasing order
  eigenVals[end] >= -1e-16 || error("negative eigenvalue(s)")

  # factor reduction
  enough = componentRetainedPercentage * sum(eigenVals)

  if componentRetainedPercentage == 1.0
    # numerical glitches might cause some factors to be discarded
    enough *= 1.1
  end

  # retain at least one factor
  components = eigenVals[1]
  retainedFactors = 1
  i = 2
  while components < enough && i <= sz
    components += eigenVals[i]
    retainedFactors += 1
    i += 1
  end

  # output is granted to have a rank <= maxrank
  retainedFactors = min(retainedFactors, maxRank)

  diagonal = zeros(sz, retainedFactors)
  for i = 1:retainedFactors
    diagonal[i, i] = sqrt(eigenVals[i])
  end

  result = jd[:vectors] * diagonal

  normalize_pseudo_root!(matrix, result)

  return result
end

function ql_norm_squared(v::Matrix{Float64}, row::Int)
  x = 0.0
  for i = 1:size(v)[2]
    x += v[row, i] * v[row, i]
  end

  return x
end

function ql_norm(v::Matrix{Float64}, row::Int)
  return sqrt(ql_norm_squared(v, row))
end

function inner_product(v::Matrix{Float64}, row1::Int, w::Matrix{Float64}, row2::Int)
  x = 0.0
  for i = 1:size(v)[2]
    x += v[row1, i] * w[row2, i]
  end

  return x
end

type OrthogonalProjection
  originalVectors::Matrix{Float64}
  multiplierCutoff::Float64
  numberVectors::Int
  numberValidVectors::Int
  dimension::Int
  validVectors::BitArray{1}
  projectedVectors::Vector{Vector{Float64}}
  orthoNormalizedVectors::Matrix{Float64}
end

function OrthogonalProjection(originalVectors::Matrix{Float64}, multiplierCutoff::Float64, tolerance::Float64)
  numberVectors, dimension = size(originalVectors)
  validVectors = trues(numberVectors)
  orthoNormalizedVectors = Matrix{Float64}(numberVectors, dimension)
  projectedVectors = Vector{Vector{Float64}}()

  currentVector = Vector{Float64}(dimension)
  for j in eachindex(validVectors)
    if validVectors[j]
      for k = 1:numberVectors
        for m = 1:dimension
          orthoNormalizedVectors[k, m] = originalVectors[k, m]
        end

        if k != j && validVectors[k]
          for l = 1:k
            if validVectors[l] && l != j
              dotProduct = inner_product(orthoNormalizedVectors, k, orthoNormalizedVectors, l)
              for n = 1:dimension
                orthoNormalizedVectors[k, n] -= dotProduct * orthoNormalizedVectors[l, n]
              end
            end
          end

          normBeforeScaling = ql_norm(orthoNormalizedVectors, k)

          if normBeforeScaling < tolerance
            validVectors[k] = false
          else
            normBeforeScalingRecip = 1.0/normBeforeScaling
            for m = 1:dimension
              orthoNormalizedVectors[k, m] *= normBeforeScalingRecip
            end
          end # end of else (norm < tolerance)
        end # end of if k != j && validVectors[k]
      end # end of for k = 1:numberVectors
      # we now have an o.n. basis for everything except j
      prevNormSquared = ql_norm_squared(originalVectors, j)

      for r = 1:numberVectors
        if validVectors[r] && r != j
          dotProduct = inner_product(orthoNormalizedVectors, j, orthoNormalizedVectors, r)
          for s = 1:dimension
            orthoNormalizedVectors[j, s] -= dotProduct * orthoNormalizedVectors[r, s]
          end
        end
      end

      projectionOnOriginalDirection = inner_product(originalVectors, j, orthoNormalizedVectors, j)
      sizeMultiplier = prevNormSquared / projectionOnOriginalDirection

      if abs(sizeMultiplier) < multiplierCutoff
        for t = 1:dimension
          currentVector[t] = orthoNormalizedVectors[j, t] * sizeMultiplier
        end
      else
        validVectors[j] = false
      end
    end # end of validVectors[j]

    push!(projectedVectors, copy(currentVector))
  end # end of j loop

  numberValidVectors = 0
  for i in eachindex(validVectors)
    numberValidVectors += validVectors[i] ? 1 : 0
  end

  return OrthogonalProjection(originalVectors, multiplierCutoff, numberVectors, numberValidVectors, dimension, validVectors, projectedVectors, orthoNormalizedVectors)
end

get_vector(op::OrthogonalProjection, i::Int) = op.projectedVectors[i]


text:
using MIToS.Utils

# Mappings
# ========

function _fill_aln_seq_ann!(aln, seq_ann::Vector{ASCIIString}, seq::Vector{UInt8},
                            init::Int, nres::Int, i)
    if length(seq) != nres
        throw(ErrorException(string("There is and aligned sequence with different number of columns [ ",
                                    length(seq), " != ", nres, " ]:\n", ascii(seq))))
    end
    @inbounds for j in 1:nres
        res = seq[j]
        aln[j,i] = Residue( res )
        if res != UInt8('-') && res != UInt8('.')
            seq_ann[j] = string(init)
            init += 1
        else
            seq_ann[j] = ""
        end
    end
    join(seq_ann, ',')
end

function _to_msa_mapping(sequences::Array{ASCIIString,1})
    nseq = size(sequences,1)
    nres = length(sequences[1])
    aln = Array(Residue,nres,nseq)
    mapp = Array(ASCIIString, nseq)
    seq_ann = Array(ASCIIString, nres)
    for i in 1:nseq
        mapp[i] = _fill_aln_seq_ann!(aln, seq_ann, sequences[i].data, 1, nres, i)
    end
    (aln', mapp)
end

function _to_msa_mapping(sequences::Array{ASCIIString,1}, ids::Array{ASCIIString,1})
    nseq = size(sequences,1)
    nres = length(sequences[1])
    aln = Array(Residue,nres,nseq)
    mapp = Array(ASCIIString, nseq)
    seq_ann = Array(ASCIIString, nres)
    sep = r"/|-"
    for i in 1:nseq
        fields = split(ids[i],sep)
        init = length(fields) == 3 ? parse(Int, fields[2]) : 1
        mapp[i] = _fill_aln_seq_ann!(aln, seq_ann, sequences[i].data, init, nres, i)
    end
    (aln', mapp)
end

# Delete Full of Gap Columns
# ==========================

"""
Deletes columns with 100% gaps, this columns are generated by inserts.
"""
function deletefullgapcolumns!(msa::AbstractMultipleSequenceAlignment, annotate::Bool=true)
    mask = columngapfraction(msa) .!= one(Float64)
    number = sum(~mask)
    if number != 0
        annotate && annotate_modification!(msa, string("deletefullgaps!  :  Deletes ", number," columns full of gaps (inserts generate full gap columns on MIToS because lowercase and dots are not allowed)"))
        filtercolumns!(msa, mask, annotate)
    end
    msa
end

function deletefullgapcolumns(msa::Matrix{Residue})
    mask = columngapfraction(msa) .!= one(Float64)
    number = sum(~mask)
    if number != 0
        return(filtercolumns(msa, mask))
    end
    msa
end

@doc """
`parse(io::Union{IO, AbstractString}, format[, output; generatemapping::Bool=false, useidcoordinates::Bool=false, deletefullgaps::Bool=true ])`

The keyword argument `generatemapping` (`false` by default) indicates if the mapping of the sequences ("SeqMap") and columns ("ColMap")
and the number of columns in the original MSA ("NCol") should be generated and saved in the annotations.
If `useidcoordinates` is `true` (default: `false`) the sequence IDs of the form "ID/start-end" are parsed and used for determining the start and end positions when the mappings are generated.
`deletefullgaps` (`true` by default) indicates if columns 100% gaps (generally inserts from a HMM) must be removed from the MSA.
By default, the ambiguous or not standard residues are replaced by gaps (i.e. `J` for leucine or isoleucine).
But, if the keyword argument `checkalphabet` is `true` (`false` by default), the sequences with residues that do not belong to the defined alphabet are deleted.""" parse

# Keepinserts
# ===========

"""
Function to keep insert columns in `parse`.
It uses the first sequence to generate the "Aligned" annotation,
and after that, convert all the characters to uppercase.
"""
function _keepinserts!(SEQS, annot)
    aligned = map(SEQS[1]) do char
        isupper(char) || char == '-' ? '1' : '0'
    end
    setannotcolumn!(annot, "Aligned", aligned)
    map!(uppercase, SEQS)
end

# Delete sequences with ambiguous or not standard residues
# ========================================================

const _alphabet = "ARNDCQEGHILKMFPSTWYV-."

function _checkalphabet(seqs::Vector{ASCIIString})
    N = length(seqs)
    use = trues(N)
    for i in 1:N
        @inbounds seq = seqs[i]
        for char in seq
            if !( uppercase(char) in _alphabet )
                @inbounds use[i] = false
                break
            end
        end
    end
    use
end

function deletenotalphabetsequences!(msa::AbstractMultipleSequenceAlignment, seqs::Vector{ASCIIString}, annotate::Bool=true)
    mask = _checkalphabet(seqs)
    number = sum(~mask)
    if number != 0
        annotate && annotate_modification!(msa, string("deletenotalphabetsequences!  :  Deletes ", number," sequences with ambiguous or not standard residues (Alphabet: ", _alphabet, " )"))
        filtersequences!(msa, mask, annotate)
    end
    msa
end

function deletenotalphabetsequences(msa::Matrix{Residue}, seqs::Vector{ASCIIString})
    mask = _checkalphabet(seqs)
    number = sum(~mask)
    if number != 0
        return(filtersequences(msa, mask))
    end
    msa
end

function _strings_to_msa(seqs::Vector{ASCIIString}, deletefullgaps::Bool, checkalphabet::Bool)
    msa = convert(Matrix{Residue}, seqs)
    if deletefullgaps && (!checkalphabet)
        return( deletefullgapcolumns(msa) )
    elseif deletefullgaps && checkalphabet
        return( deletefullgapcolumns( deletenotalphabetsequences(msa, seqs) ) )
    elseif (!deletefullgaps) && checkalphabet
        return( deletenotalphabetsequences(msa, seqs) )
    else
        return(msa)
    end
end


text:
#physics sims starting with incline planes

using Plots

struct Plane
	muS::Float64
	muK::Float64
	theta::Float64
end

function makeDefault(p::Plane)
	plt = plot(xlims=[0,100],ylims=[0,100],legend=false)
	plot!([100-100/(tand(p.theta)),100],[100,0])
	plt
end

function makeDefault(ps::Array{Plane,1})
	plt = plot(xlims=[0,100],ylims=[0,100],legend=false)
	for p in ps
		plot!([100-100/(tand(p.theta)),100],[100,0])
	end
	plt
end

function addBox(x0,p::Plane,title::String)
	plt=makeDefault(p)
	box = [((100-x0)*cosd(p.theta),x0*sind(p.theta)),((100-(x0+10))*cosd(p.theta),(x0+10)*sind(p.theta))]
	plot!(box,w=10,title="μK=$(p.muK), μS=$(p.muS), θ=$(p.theta) $(title)")
	plt
end

#boxes at the same distance x0 from the end
function addBox(x0::Array{Float64},ps::Array{Plane,1},title::String)
	plt=makeDefault(ps)
	i=1
	for p in ps
		box = [(100-x0[i]*cosd(p.theta),x0[i]*sind(p.theta)),(100-(x0[i]+10)*cosd(p.theta),(x0[i]+10)*sind(p.theta))]
		plot!(box,w=10)
		i+=1;
	end
	title!(title)
	plt
end

function travelPlane(x0,ps::Array{Plane,1})
	plt = makeDefault(ps)
	thetas = sort!([p.theta for p in ps];rev=true)
	t_end = sqrt(3*abs(x0)*sind(thetas[1])/9.8)
	vals = range(0.0,t_end,length = 100)
	anim = @animate for t in vals
		for p in ps
			y=x0-4.9*sin(p.theta)*t^2
			box = [(100-y*cosd(p.theta),y*sind(p.theta)),(100-(y+10)*cosd(p.theta),(y+10)*sind(p.theta))]
			plot!(box,w=10)
		end

	end
	gif(anim,"boxOnPlanes.gif",fps=25)
end
function travelPlane1(x0,ps::Array{Plane,1})
	plt = makeDefault(ps)
	thetas = sort!([p.theta for p in ps];rev=true)
	t_end = sqrt(2*abs(x0)*sind(thetas[1])/9.8)
	vals = range(0.0,t_end,length = 100)
	anim = @animate for t in vals
		addBox([x0-(4.9*t^2)*sind(p.theta) for p in ps],ps,"t=$t seconds")
	end
	gif(anim,"boxOnPlanes.gif",fps=25)
end

function travelPlane(x0,p::Plane)

	t_end = sqrt(3*abs(x0)*sind(p.theta)/9.8)

	vals = range(0.0,t_end,length = 100)
	anim = @animate for t in vals
		addBox(x0-4.9*sin(p.theta)*t^2,p,"t=$(t) seconds")
	end
	gif(anim,"boxPlane.gif",fps=25)
end

max_stars_repo_name: thomastjdavis/PhysicsDemos

max_issues_repo_name: thomastjdavis/PhysicsDemos

max_forks_repo_name: thomastjdavis/PhysicsDemos


text:
@testset "OLS" begin
    T, N, K = 100, 2, 4
    X = randn(T, K)
    b = randn(K, N)
    Y = X * b
    m = OLS(Y, X)
    @test m.coef ≈ b
    @test modelmatrix(m) === X
    @test coef(m) === m.coef
    @test residuals(m) === m.resid
    @test sprint(show, m) == "OLS regression"
end

@testset "LeastSquaresLP" begin
    e = LeastSquaresLP()
    @test sprint(show, e) == "LeastSquaresLP"
    @test sprint(show, MIME("text/plain"), e) == "Ordinary least squares local projection"
end

@testset "LeastSquaresLPResult" begin
    er = LeastSquaresLPResult([OLS(rand(5), rand(5,2))])
    @test sprint(show, er) == "LeastSquaresLPResult"
end

@testset "LocalProjectionResult" begin
    B = cat(randn(5,2,1), randn(5,2,1); dims=3)
    V = cat(randn(10,10), randn(10,10); dims=3)
    r = LocalProjectionResult(B, V, [1,2], LeastSquaresLP(), nothing, HRVCE(), VarName[:y1, :y2], VarName[:x], VarName[:w1, :w2], Dict{VarName,Int}(:y1=>1, :y2=>2), Dict{VarName,Int}(:x=>1), Dict{VarName,Int}(:w1=>1,:w2=>2,:y1=>3,:y2=>4), 2, 0, nothing, nothing, nothing, nothing, nothing, nothing, false, true)
    @test coef(r, 1, :x) == B[1,1,1]
    @test coef(r, 2, :w1, yname=:y2, lag=1) == B[2,2,2]
    @test coef(r, 2, :w1, lag=2) == B[4,1,2]
    @test vcov(r, 1, :x) == V[1,1,1]
    @test vcov(r, 2, :x, yname1=:y2, xwname2=:w1, lag2=1) == V[2,4,2]
    @test vcov(r, 2, :x, yname1=:y2, xwname2=:w1, yname2=:y1, lag2=1) == V[2,3,2]
    @test vcov(r, 2, :w1, yname1=:y2, lag1=2, xwname2=:w1, yname2=:y1, lag2=1) == V[8,3,2]
    @test vcov(r, 2, :w2, yname1=:y2, lag1=2, xwname2=:w1, yname2=:y1) == V[10,7,2]

    @test sprint(show, r) == "LocalProjectionResult"
    @test sprint(show, MIME("text/plain"), r) == """
        LocalProjectionResult with 2 lags over 2 horizons:
        ──────────────────────────────────────────────────────────────────────────────
        Variable Specifications
        ──────────────────────────────────────────────────────────────────────────────
        Outcome variables:              y1 y2    Minimum horizon:                    0
        Regressor:                          x    Lagged controls:                w1 w2
        ──────────────────────────────────────────────────────────────────────────────"""
    r = LocalProjectionResult(ones(1,1,1), ones(1,1,1), [1], LeastSquaresLP(), nothing, HRVCE(), VarName[:y], VarName[], VarName[:w], Dict{VarName,Int}(:y=>1), Dict{VarName,Int}(), Dict{VarName,Int}(:w=>1), 1, 0, nothing, nothing, nothing, nothing, nothing, nothing, false, true)
    @test sprint(show, MIME("text/plain"), r) == """
        LocalProjectionResult with 1 lag over 1 horizon:
        ──────────────────────────────────────────────────────────────────────────────
        Variable Specifications
        ──────────────────────────────────────────────────────────────────────────────
        Outcome variable:                   y    Minimum horizon:                    0
        Regressor:                               Lagged control:                     w
        ──────────────────────────────────────────────────────────────────────────────"""
    r = LocalProjectionResult(ones(1,1,1), ones(1,1,1), [1], LeastSquaresLP(), nothing, HRVCE(), VarName[Cum(:y)], VarName[Cum(:x)], VarName[:w], Dict{VarName,Int}(Cum(:y)=>1), Dict{VarName,Int}(Cum(:x)=>1), Dict{VarName,Int}(:w=>1), 1, 0, nothing, nothing, nothing, nothing, VarName[Cum(:x)], VarName[:z1,:z2], false, true)
    @test sprint(show, MIME("text/plain"), r) == """
        LocalProjectionResult with 1 lag over 1 horizon:
        ──────────────────────────────────────────────────────────────────────────────
        Variable Specifications
        ──────────────────────────────────────────────────────────────────────────────
        Outcome variable:              Cum(y)    Minimum horizon:                    0
        Regressor:                     Cum(x)    Lagged control:                     w
        Endogenous variable:           Cum(x)    Instruments:                    z1 z2
        ──────────────────────────────────────────────────────────────────────────────"""
end

@testset "_makeYX" begin
    T = 100
    ys = [randn(T), randn(T)]
    xs = [randn(T), randn(T)]
    ws = ys
    Y, X, T1, eT = _makeYX(ys, xs, ws, nothing, 3, 5, nothing)
    @test T1 == T-8
    @test size(Y) == (T1, 2)
    @test size(X) == (T1, 8)
    @test Y[:,1] == ys[1][9:end]
    @test Y[:,2] == ys[2][9:end]
    @test X[:,1] == xs[1][4:end-5]
    @test X[:,3] == ws[1][3:end-6]
    @test X[:,4] == ws[2][3:end-6]
    @test X[:,8] == ws[2][1:end-8]
    @test all(eT)

    sts = (rand(T), rand(T))
    Y, X, T1, eT = _makeYX(ys, xs, ws, sts, 1, 0, nothing)
    @test size(X) == (99, 6)
    @test X[:,1] == xs[1][2:end,1]
    @test X[:,2] == xs[2][2:end,1]
    @test X[:,3] == ys[1][1:99,1].*sts[1][2:end]
    @test X[:,4] == ys[1][1:99,1].*sts[2][2:end]
    @test X[:,5] == ys[2][1:99,1].*sts[1][2:end]
    @test X[:,6] == ys[2][1:99,1].*sts[2][2:end]
    @test all(eT)

    ys = (randn(T),)
    xs = ()
    ws = ys
    Y, X, T1, eT = _makeYX(ys, xs, ws, nothing, 1, 0, nothing)
    @test T1 == 99
    @test size(Y) == (T1, 1)
    @test size(X) == (T1, 1)
    @test Y == reshape(ys[1][2:end], T1, 1)
    @test X == reshape(ws[1][1:end-1], T1, 1)
    @test all(eT)

    ys[1][2], ys[1][3] = NaN, Inf
    xs = (convert(Vector{Union{Float64, Missing}}, randn(T)),)
    xs[1][3] = missing
    Y, X, T1, eT = _makeYX(ys, xs, ws, nothing, 1, 0, nothing)
    @test T1 == 96
    @test size(Y) == (T1, 1)
    @test size(X) == (T1, 2)
    @test Y == reshape(ys[1][5:end], T1, 1)
    @test X[:,1] == xs[1][5:end]
    @test X[:,2] == ws[1][4:end-1]
    @test eT == ((1:99).>3)

    Y, X, T1, eT = _makeYX(ys, xs, ws, nothing, 1, 1, (1:100).<=90)
    @test T1 == 85
    @test size(Y) == (T1, 1)
    @test size(X) == (T1, 2)
    @test Y == reshape(ys[1][6:90], T1, 1)
    @test X[:,1] == xs[1][5:89]
    @test X[:,2] == ws[1][4:88]
    # eT covers 2:99 in full data
    @test eT == ((1:98).>3) .& ((1:98).<=88)

    Y, X, T1, eT = _makeYX(ys, xs, ws, nothing, 1, 1, ((1:100).<60).|((1:100).>=70))
    @test T1 == 83
    @test Y[:] == vcat(ys[1][6:59], ys[1][72:end])
    @test X[:,1] == vcat(xs[1][5:58], xs[1][71:end-1])
    @test X[:,2] == vcat(ws[1][4:57], ws[1][70:end-2])
    @test eT == ((1:98).>3) .& (((1:98).<58) .| ((1:98).>=70))

    ys = ()
    @test_throws ArgumentError _makeYX(ys, xs, ws, nothing, 1, 0, nothing)
    ws = ()
    @test_throws ArgumentError _makeYX(ys, xs, ws, nothing, 1, 0, nothing)
    @test_throws ArgumentError _makeYX(ys, xs, ws, nothing, 0, 0, nothing)

    ys = ([Inf, NaN, 1.0],)
    xs = ()
    ws = ys
    @test_throws ArgumentError _makeYX(ys, xs, ws, nothing, 1, 0, nothing)
end

@testset "lp" begin
    T = 100
    ys = [randn(T), randn(T)]
    xs = [randn(T), randn(T)]
    ws = ys
    b, v, t, m = _lp(ys, xs, ws, nothing, 3, 5, SimpleVCE(), nothing)
    @test size(b) == (8, 2)
    @test size(v) == (16, 16)
    @test t == 92

    # Reproduce estimates from example based on Gertler & Karadi (2015)
    # Specification is non-augmented local projection
    # EHW standard errors
    # https://github.com/jm4474/Lag-augmented_LocalProjections/blob/master/examples/gk.m
    # The intercept estimated will be the same if not detrend the level before estimation
    # The left division in Matlab code uses QR factorization to solve OLS
    # The use of Cholesky factorization results in some trivial discrepancy
    df = exampledata(:gk)
    ns = (:logcpi, :logip, :ff, :ebp, :ff4_tc)
    r = lp(df, :ebp, wnames=ns, nlag=12, nhorz=48, vce=HRVCE())
    @test coef(r, 1, :ff4_tc, lag=1) ≈ 0.394494467439933 atol=1e-8
    @test coef(r, 5, :ff4_tc, lag=1) ≈ 0.147288229090455 atol=1e-8
    @test sqrt(vcov(r, 1, :ff4_tc, lag1=1)) ≈ 0.473834526724786 atol=1e-8
    @test sqrt(vcov(r, 5, :ff4_tc, lag1=1)) ≈ 0.737092370022356 atol=1e-8
    # Compare additional estimates obtained from Matlab lp before being turned into irf
    # Run in gk.m: lp(Y, 11, 1, 4, false, false)
    @test coef(r, 1, :constant) ≈ -2.945365821418727 atol=1e-8
    @test coef(r, 1, :logip, lag=1) ≈ -0.060094199554036 atol=1e-8
    @test vcov(r, 1, :constant) ≈ 0.641933721888132 atol=1e-8
    @test vcov(r, 1, :logcpi, lag1=1) ≈ 0.007331238288153 atol=1e-8
    @test vcov(r, 1, :logcpi, lag1=12) ≈ 0.004780116443880 atol=1e-8

    f = irf(r, :ebp, :ff4_tc, lag=1)
    @test coef(f)[1] ≈ 0.394494467439933 atol=1e-8
    @test coef(f)[20] ≈ -0.552849402428304 atol=1e-8
    @test stderror(f)[1] ≈ 0.473834526724786 atol=1e-8
    @test stderror(f)[20] ≈ 1.005046789996172 atol=1e-8
    ci = confint(f)
    @test ci[1][1] ≈ -0.384893972486403 atol=1e-8
    @test ci[2][1] ≈ 1.173882907366730 atol=1e-8
    @test ci[1][20] ≈ -2.206004261810359 atol=1e-8
    @test ci[2][20] ≈ 1.100305456953579 atol=1e-8

    ns1 = [3, :logip, 5, :ff4_tc]
    r1 = lp(df, 11, wnames=ns1, nlag=12, nhorz=1, vce=HRVCE())
    @test r1.B[1:3,1,1] ≈ r.B[1:3,1,1]
    @test r1.V[1:3,1:3,1] ≈ r.V[1:3,1:3,1]

    rn = lp(df, :ebp, xnames=:ff4_tc, wnames=(:ff4_tc,), nlag=12, nhorz=2,
        normalize=:ff4_tc=>:ff, vce=HRVCE())
    riv = lp(df, :ebp, xnames=:ff, wnames=(:ff4_tc,), nlag=12, nhorz=2,
        iv=:ff=>:ff4_tc, vce=HRVCE())
    @test coef(rn, 1, :ff4_tc) ≈ coef(riv, 1, :ff)
    @test coef(rn, 2, :ff4_tc) ≈ coef(riv, 2, :ff)
    @test vcov(rn, 1, :ff4_tc) ≈ vcov(riv, 1, :ff)
    @test vcov(rn, 2, :ff4_tc) ≈ vcov(riv, 2, :ff)
    @test rn.normnames == [:ff4_tc]
    @test rn.normtars == [:ff]
    @test rn.normmults ≈ [-5.3028242533068495] atol=1e-8
    @test riv.endonames == VarName[:ff]
    @test riv.ivnames == VarName[:ff4_tc]

    @test_throws ArgumentError lp(df, :ebp, xnames=:ff4_tc, normalize=:no=>:ff)
    @test_throws ArgumentError lp(df, :ebp, xnames=:ff, iv=:no=>:ff4_tc)
    @test_throws ArgumentError lp(df, :ebp, xnames=:ff, iv=:ff=>())
    @test_throws ArgumentError lp(df, :ebp, xnames=:ff, iv=:ff=>1.0)

    # Compare standard errors based on EWC with Matlab results
    r2 = lp(r, HARVCE(EWC()))
    f2 = irf(r2, :ebp, :ff4_tc, lag=1)
    @test f2.B == f.B
    @test f2.T == f.T
    @test stderror(f2)[1] ≈ 0.408141982124863 atol=1e-8
    @test stderror(f2)[20] ≈ 1.037575254556690 atol=1e-8
    ci2 = confint(f2)
    @test ci2[1][1] ≈ -0.318073956754668 atol=1e-8
    @test ci2[2][1] ≈ 1.107062891634995 atol=1e-8
    @test ci2[1][20] ≈ -2.371771071486807 atol=1e-8
    @test ci2[2][20] ≈ 1.266072266630027 atol=1e-8

    # Reproduce point estimates in Table 1 from Ramey and Zubairy (2018)
    # Compare estimates produced by ivreg2 in Stata replication files
    # `jordagk.do` and `jordagk_twoinstruments.do`
    # Need to run `set type double` before executing the Stata do files for more precision
    # Use `bro h multlin1` to see the Stata results
    df = exampledata(:rz)

    # Baseline specifications
    r1 = lp(df, Cum(:y), xnames=Cum(:g), wnames=(:newsy, :y, :g), iv=Cum(:g)=>:newsy,
        nlag=4, nhorz=17, addylag=false, firststagebyhorz=true)
    f1 = irf(r1, Cum(:y), Cum(:g))
    @test coef(f1)[1] ≈ 1.306459420753 atol=1e-8
    @test coef(f1)[9] ≈ .6689611146745 atol=1e-9
    @test coef(f1)[17] ≈ .7096110638625 atol=1e-9
    # Standard errors based on EWC rather than Newey-West in original paper
    @test stderror(f1)[1] ≈ 0.4286854811143802 atol=1e-8
    @test stderror(f1)[9] ≈ 0.10921291172237432 atol=1e-8
    @test stderror(f1)[17] ≈ 0.20049776900175723 atol=1e-8

    r2 = lp(df, Cum(:y), xnames=Cum(:g), wnames=(:newsy, :y, :g), iv=Cum(:g)=>(:newsy, :g),
        nlag=4, nhorz=16, minhorz=1, addylag=false, firststagebyhorz=true)
    f2 = irf(r2, Cum(:y), Cum(:g))
    @test coef(f2)[1] ≈ .218045283661 atol=1e-9
    @test coef(f2)[8] ≈ .4509452247638 atol=1e-9
    @test coef(f2)[16] ≈ .5591002357465 atol=1e-9
    # Standard errors based on EWC rather than Newey-West in original paper
    @test stderror(f2)[1] ≈ 0.16423131478595118 atol=1e-8
    @test stderror(f2)[8] ≈ 0.10754779358937377 atol=1e-8
    @test stderror(f2)[16] ≈ 0.14263798115850218 atol=1e-8

    # Omit WWII
    r1 = lp(df, Cum(:y), xnames=Cum(:g), wnames=(:newsy, :y, :g), iv=Cum(:g)=>:newsy,
        nlag=4, nhorz=17, addylag=false, firststagebyhorz=true, subset=df.wwii.==0)
    f1 = irf(r1, Cum(:y), Cum(:g))
    @test coef(f1)[1] ≈ 1.7188411171 atol=1e-8
    @test coef(f1)[9] ≈ .7672214331391 atol=1e-9
    @test coef(f1)[17] ≈ .7167179685971 atol=1e-9

    r2 = lp(df, Cum(:y), xnames=Cum(:g), wnames=(:newsy, :y, :g), iv=Cum(:g)=>(:newsy, :g),
        nlag=4, nhorz=16, minhorz=1, addylag=false, firststagebyhorz=true, subset=df.wwii.==0)
    f2 = irf(r2, Cum(:y), Cum(:g))
    @test coef(f2)[1] ≈ -.0342369278527 atol=1e-10
    @test coef(f2)[8] ≈ .2570200212481 atol=1e-9
    @test coef(f2)[16] ≈ .2225592580744 atol=1e-9

    # State dependency
    # Compare estimates generated from `jordagk_twoinstruments.do`
    # Use `bro multexpbh* multrecbh*` to see the Stata results

    # Construct the lag of slack
    df.rec = similar(df.slack)
    df.rec[1] = missing
    # Fill the beginning missing values
    df.rec[2:5] .= true
    df.rec[6:end] .= view(df.slack, 5:length(df.rec)-1)
    df.exp = 1.0 .- df.rec
    df.recnewsy = df.rec.*df.newsy
    df.expnewsy = df.exp.*df.newsy
    df.recg = df.rec.*df.g
    df.expg = df.exp.*df.g

    # nomit
    r2 = lp(df, Cum(:y), xnames=(Cum(:g,:rec), Cum(:g,:exp), :rec), wnames=(:newsy, :y, :g),
        iv=(Cum(:g,:rec), Cum(:g,:exp))=>(:recnewsy, :expnewsy, :recg, :expg),
        states=(:rec, :exp), nlag=4, nhorz=16, minhorz=1, addylag=false, firststagebyhorz=true)
    f2rec = irf(r2, Cum(:y), Cum(:g,:rec))
    @test coef(f2rec)[1] ≈ .2718093668839 atol=1e-9
    @test coef(f2rec)[8] ≈ .6357587189621 atol=1e-9
    @test coef(f2rec)[16] ≈ .6784987089259 atol=1e-9
    f2exp = irf(r2, Cum(:y), Cum(:g,:exp))
    @test coef(f2exp)[1] ≈ .2660710337235 atol=1e-9
    @test coef(f2exp)[8] ≈ .3511868388287 atol=1e-9
    @test coef(f2exp)[16] ≈ .373442191223 atol=1e-9

    # wwii
    r2 = lp(df, Cum(:y), xnames=(Cum(:g,:rec), Cum(:g,:exp), :rec), wnames=(:newsy, :y, :g),
        iv=(Cum(:g,:rec), Cum(:g,:exp))=>(:recnewsy, :expnewsy, :recg, :expg),
        states=(:rec, :exp), nlag=4, nhorz=16, minhorz=1, addylag=false, firststagebyhorz=true,
        subset=df.wwii.==0)
    f2rec = irf(r2, Cum(:y), Cum(:g,:rec))
    @test coef(f2rec)[1] ≈ .345617820757 atol=1e-9
    @test coef(f2rec)[8] ≈ 1.350903038984 atol=1e-8
    @test coef(f2rec)[16] ≈ 1.365148452565 atol=1e-8
    f2exp = irf(r2, Cum(:y), Cum(:g,:exp))
    @test coef(f2exp)[1] ≈ .010365014606 atol=1e-10
    @test coef(f2exp)[8] ≈ .2171636455168 atol=1e-9
    @test coef(f2exp)[16] ≈ .2140385336746 atol=1e-9

    @test_logs (:warn, "firststagebyhorz=false while endogenous variables contain Cum")
        lp(df, Cum(:y), xnames=Cum(:g), wnames=(:y,), iv=Cum(:g)=>:newsy, addylag=false)
    @test_logs (:warn, "addylag=true while outcome variables contain Cum")
        lp(df, Cum(:y), xnames=Cum(:g))
end


text:
# Shrinking functions for default data types

# shrinks `true` to `Bool[false]`` and `false` to `Bool[]`
function shrink(t::Bool)
    if t
        return Bool[false]
    else
        return Bool[]
    end
end

function shrink(w::T) where T <: Unsigned
    ret = T[ w & ~(one(T) << mask) for mask in zero(T):(sizeof(T)*8 - 1) ]
    w > zero(T) && pushfirst!(ret, w - 0x1)
    w > one(T) && pushfirst!(ret, w - 0x2)
    pushfirst!(ret, w >> 0x1)
    return sort!(filter!(!=(w), unique!(ret)))
end

# shrinks a signed value by shrinking its absolute value like an unsigned and then negating that
function shrink(w::T) where T <: Signed
    ret = signed.(shrink(unsigned(abs(w))))
    append!(ret, ret)
    retView = @view ret[end÷2+1:end]
    map!(-, retView, retView)
    return filter!(unique!(ret)) do x
        x != w && x != -w
    end
end

# shrinks a character by shrinking its codepoint
function shrink(w::Char)
    c = codepoint(w)
    ret = Char[]
    c > 2 && push!(ret, c ÷ 2)
    c > 0 && push!(ret, c - 1)
    return ret
end

# drops a character and shrinks a character to form new strings
function shrink(s::String)
    ret = String[]
    io = IOBuffer()
    for i in eachindex(s)
        head = @view s[begin:i-1]
        tail = @view s[i+1:end]
        write(io, head)
        write(io, tail)
        push!(ret, String(take!(io)))

        shrinks = shrink(s[i])
        for sh in shrinks
            write(io, head)
            write(io, sh)
            write(io, tail)
            push!(ret, String(take!(io)))
        end
    end
    return ret
end

###########
# Vectors & Tuples
###########

# drop each index once
drops(root) = (deleteat!(deepcopy(root), i) for i in eachindex(root) if length(root) > 1)

# create all tuples with all identity, except for one place which we shrink
allFuncs(root) = (( i == idx ? shrink : (x -> [x]) for i in 1:length(root) ) for idx in eachindex(root))

# apply the shrinking functions and filter unsuccessful shrinks
shrunkEls(root) = ifilter((x -> !(any(isempty, x) || all(==(root), x))),
                    ((f(x) for (f,x) in zip(funcs, root)) for funcs in allFuncs(root)))

getProd(a, ::Type{<:Tuple}) = a
getProd(a, ::Type{<:Vector}) = [a...]
# we have a number of shrink results, take the product of all for all shrink combinations
prods(root::T) where T = flatten((getProd(p, T) for p in iproduct(shrunks...)) for shrunks in shrunkEls(root))

# finally, filter out productions that resulted back in the original root for efficiency
shrinks(root) = (ifilter(!=(root), prods(root)))

# combine with drops and we got lazy shrinks (～￣▽￣)～
shrink(rootEl::Vector) = iunique(flatten((drops(rootEl), shrinks(rootEl))))
# and we can reuse the methods for tuples as well
shrink(rootEl::Tuple) = iunique(shrinks(rootEl))


text:
@testset "Accelerated Gradient Descent" begin
    f(x) = x[1]^4
    function g!(storage, x)
        storage[1] = 4 * x[1]^3
        return
    end

    initial_x = [1.0]
    options = Optim.Options(show_trace = debug_printing, allow_f_increases=true)
    results = Optim.optimize(f, g!, initial_x, AcceleratedGradientDescent(), options)
    @test norm(Optim.minimum(results)) < 1e-6
    @test summary(results) == "Accelerated Gradient Descent"

    run_optim_tests(AcceleratedGradientDescent(); skip = ("Large Polynomial","Parabola"),
                                                  convergence_exceptions = (("Rosenbrock", 1),("Rosenbrock", 2)),
                                                  iteration_exceptions = (("Powell", 1100),
                                                                          ("Polynomial", 1500)))
end


text:
## wiener filter

@doc doc"""
    argmaxf_lnP(ϕ,                ds::DataSet; kwargs...)
    argmaxf_lnP(ϕ, θ, ds::DataSet; kwargs...)
    argmaxf_lnP(Lϕ,               ds::DataSet; kwargs...)

Computes either the Wiener filter at fixed $\phi$, or a sample from this slice
along the posterior.

Keyword arguments:

* `which` — `:wf`, `:sample`, or `fluctuation` to compute 1) the Wiener filter,
  i.e. the best-fit of $\mathcal{P}(f\,|\,\phi,d)$, 2) a sample from
  $\mathcal{P}(f\,|\,\phi,d)$, or 3) a sample minus the Wiener filter, i.e. the
  fluctuation on top of the mean.
* `fstart` — starting guess for `f` for the conjugate gradient solver
* `conjgrad_kwargs` — Passed to the inner call to [`conjugate_gradient`](@ref)

"""
argmaxf_lnP(ϕ::Field,                ds::DataSet; kwargs...) = argmaxf_lnP(cache(ds.L(ϕ),ds.d), NamedTuple(), ds; kwargs...)
argmaxf_lnP(ϕ::Field, θ, ds::DataSet; kwargs...) = argmaxf_lnP(cache(ds.L(ϕ),ds.d), θ,            ds; kwargs...)

function argmaxf_lnP(
    Lϕ,
    θ,
    ds::DataSet;
    which = :wf,
    fstart = nothing,
    preconditioner = :diag,
    conjgrad_kwargs = (tol=1e-1,nsteps=500)
)

    @unpack d, Cn, Cn̂, Cf, M, M̂, B, B̂ = ds(θ)

    Δ = d - nonCMB_data_components(θ,ds)
    b = 0
    if (which in (:wf, :sample))
        b += Lϕ'*B'*M'*(Cn\Δ)
    end
    if (which in (:fluctuation, :sample))
        ξf = simulate(Cf; d.Nbatch)
        ξn = simulate(Cn; d.Nbatch)
        b += Cf\ξf + Lϕ'*B'*M'*(Cn\ξn)
    end

    A_diag  = pinv(Cf) +     B̂'*M̂'*pinv(Cn̂)*M̂*B̂
    A_zeroϕ = pinv(Cf) +     B'*M'*pinv(Cn̂)*M*B
    A       = pinv(Cf) + Lϕ'*B'*M'*pinv(Cn)*M*B*Lϕ

    A_preconditioner = @match preconditioner begin
        :diag  => A_diag
        :zeroϕ => FuncOp(op⁻¹ = (b -> (conjugate_gradient(A_diag, A_zeroϕ, b, zero(b); conjgrad_kwargs.tol))))
        _      => error("Unrecognized preconditioner='$preconditioner'")
    end

    conjugate_gradient(A_preconditioner, A, b, (isnothing(fstart) ? zero(b) : fstart); conjgrad_kwargs...)

end

@doc doc"""
    Σ(ϕ::Field,  ds; [conjgrad_kwargs])
    Σ(Lϕ,        ds; [conjgrad_kwargs])

An operator for the data covariance, `Cn + M*B*L*Cf*L'*B'*M'`, which can
applied and inverted. `conjgrad_kwargs` are passed to the underlying call to
`conjugate_gradient`.
"""
Σ(ϕ::Field, ds; kwargs...) = Σ(ds.L(ϕ), ds; kwargs...)
function Σ(Lϕ, ds; conjgrad_kwargs=(tol=1e-1,nsteps=500))
    @unpack d,M,B,Cn,Cf,Cn̂,B̂,M̂ = ds
    SymmetricFuncOp(
        op   = x -> (Cn + M*B*Lϕ*Cf*Lϕ'*B'*M')*x,
        op⁻¹ = x -> conjugate_gradient((Cn̂ .+ M̂*B̂*Cf*B̂'*M̂'), Σ(Lϕ, ds), x; conjgrad_kwargs...)
    )
end

@doc doc"""

    MAP_joint(ds::DataSet; kwargs...)

Compute the maximum a posteriori (i.e. "MAP") estimate of the joint
posterior, $\mathcal{P}(f,\phi,\theta\,|\,d)$, or compute a
quasi-sample.

Keyword arguments:

* `nsteps` — The maximum number of iterations for the maximizer.
* `ϕstart = 0` — Starting point of the maximizer.
* `ϕtol = nothing` — If given, stop when `ϕ` updates reach this
  tolerance. `ϕtol` is roughly the relative per-pixel standard
  deviation between changes to `ϕ` and draws from the `ϕ` prior.
  Values in the range $10^{-2}-10^{-4}$ are reasonable.
* `lbfgs_rank = 5` — The maximum rank of the LBFGS approximation to
   the Hessian.
* `conjgrad_kwargs = (;)` — Passed to the inner call to
  [`conjugate_gradient`](@ref).
* `progress = true` — Whether to show the progress bar.
* `Nϕ = :qe` — Noise to use in the initial approximation to the
   Hessian. Can give `:qe` to use the quadratic estimate noise.
* `quasi_sample = false` — `false` to compute the MAP, `true` to
   iterate quasi-samples, or an integer to compute a fixed-seed
   quasi-sample.
* `history_keys` — What quantities to include in the returned
  `history`. Can be any subset of `(:f, :f°, :ϕ, :∇ϕ_lnP, :χ², :lnP)`.

Returns a tuple `(f, ϕ, history)` where `f` is the best-fit (or
quasi-sample) field, `ϕ` is the lensing potential, and `history`
contains the history of steps during the run.

"""
MAP_joint(ds::DataSet; kwargs...) = MAP_joint(NamedTuple(), ds; kwargs...)
function MAP_joint(
    θ :: NamedTuple,
    ds :: DataSet;
    nsteps = 20,
    lbfgs_rank = 5,
    Nϕ = :qe,
    ϕstart = nothing,
    fstart = nothing,
    ϕtol = nothing,
    progress::Bool = true,
    verbosity = (0,0),
    linesearch = OptimKit.HagerZhangLineSearch(verbosity=verbosity[2], maxiter=5),
    conjgrad_kwargs = (tol=1e-1,nsteps=500),
    quasi_sample = false,
    preconditioner = :diag,
    history_keys = (:lnP,),
    aggressive_gc = fieldinfo(ds.d).Nx>=1024 & fieldinfo(ds.d).Ny>=1024,
)

    dsθ = copy(ds(θ))
    dsθ.G = 1 # MAP estimate is invariant to G so avoid wasted computation

    ϕ = Map(isnothing(ϕstart) ? zero(diag(ds.Cϕ)) : ϕstart)

    # compute approximate inverse ϕ Hessian used in gradient descent, possibly
    # from quadratic estimate
    if (Nϕ == :qe)
        Nϕ = quadratic_estimate(dsθ).Nϕ/2
    end
    Hϕ⁻¹ = (Nϕ == nothing) ? dsθ.Cϕ : pinv(pinv(dsθ.Cϕ) + pinv(Nϕ))

    history = []
    argmaxf_lnP_kwargs = (
        which = (quasi_sample==false) ? :wf : :sample,
        preconditioner = preconditioner,
        conjgrad_kwargs = (history_keys=(:i,:res), progress=false, conjgrad_kwargs...),
    )
    pbar = Progress(nsteps, (progress ? 0 : Inf), "MAP_joint: ")
    ProgressMeter.update!(pbar)

    (f, argmaxf_lnP_history) = argmaxf_lnP(
        (isnothing(ϕstart) ? 1 : ϕ), θ, dsθ;
        fstart, argmaxf_lnP_kwargs...
    )
    f°, = mix(f, ϕ, dsθ)
    lastϕ = nothing
    push!(history, select((;f,f°,ϕ,∇ϕ_lnP=nothing,χ²=nothing,lnP=nothing,argmaxf_lnP_history), history_keys))

    # objective function (with gradient) to maximize
    @⌛ function objective(ϕ; need_gradient=true)
        (
            @⌛(sum(unbatch(-2lnP(:mix,f°,ϕ,dsθ)))),
            need_gradient ? @⌛(gradient(ϕ->-2lnP(:mix,f°,ϕ,dsθ), ϕ)[1]) : nothing
        )
    end
    # function to compute after each optimization iteration, which
    # recomputes the best-fit f given the current ϕ
    @⌛ function finalize!(ϕ,χ²,∇ϕ_lnP,i)
        if isa(quasi_sample,Int)
            seed!(global_rng_for(f),quasi_sample)
        end
        (f, argmaxf_lnP_history) = @⌛ argmaxf_lnP(
            ϕ, θ, dsθ;
            fstart = f,
            argmaxf_lnP_kwargs...
        )
        aggressive_gc && cuda_gc()
        f°, = mix(f, ϕ, dsθ)
        ∇ϕ_lnP .= @⌛ gradient(ϕ->-2lnP(:mix,f°,ϕ,dsθ), ϕ)[1]
        χ²s = @⌛ -2lnP(:mix,f°,ϕ,dsθ)
        χ² = sum(unbatch(χ²s))
        values = [("step",i), ("χ²",χ²s), ("Ncg",length(argmaxf_lnP_history))]
        hasproperty(linesearch, :α) && push!(values, ("α", linesearch.α))
        next!(pbar, showvalues=values)
        push!(history, select((;f,f°,ϕ,∇ϕ_lnP,χ²,lnP=-χ²/2,argmaxf_lnP_history), history_keys))
        if (
            !isnothing(ϕtol) &&
            !isnothing(lastϕ) &&
            sum(unbatch(norm(LowPass(1000) * (sqrt(ds.Cϕ) \ (ϕ - lastϕ))) / sqrt(2length(ϕ)))) < ϕtol
        )
            ∇ϕ_lnP = zero(∇ϕ_lnP) # this stops the solver here
        else
            lastϕ = ϕ
        end
        ϕ, χ², ∇ϕ_lnP
    end

    # run optimization
    ϕ, = @⌛ optimize(
        objective,
        Map(ϕ),
        OptimKit.LBFGS(
            lbfgs_rank;
            maxiter = nsteps,
            verbosity = verbosity[1],
            linesearch = linesearch
        );
        finalize!,
        inner = (_,ξ1,ξ2)->sum(unbatch(dot(ξ1,ξ2))),
        precondition = (_,η)->Map(Hϕ⁻¹*η),
    )

    ProgressMeter.finish!(pbar)

    f, ϕ, history

end

OptimKit._scale!(η::Field, β) = η .*= β
OptimKit._add!(η::Field, ξ::Field, β) = η .+= β .* ξ

@doc doc"""

    MAP_marg(ds; kwargs...)

Compute the maximum a posteriori (i.e. "MAP") estimate of the marginal posterior,
$\mathcal{P}(\phi,\theta\,|\,d)$.

"""
MAP_marg(ds::DataSet; kwargs...) = MAP_marg(NamedTuple(), ds; kwargs...)
function MAP_marg(
    θ,
    ds :: DataSet;
    ϕstart = nothing,
    Nϕ = :qe,
    nsteps = 10,
    nsteps_with_meanfield_update = 4,
    conjgrad_kwargs = (tol=1e-1,nsteps=500),
    α = 0.2,
    weights = :unlensed,
    Nsims = 50,
    Nbatch = 1,
    progress::Bool = true,
    aggressive_gc = fieldinfo(ds.d).Nx >=512 & fieldinfo(ds.d).Ny >=512
)

    mod(Nsims+1,nworkers())==0 || @warn "MAP_marg is most efficient when Nsims+1 is divisible by the number of workers." maxlog=1

    ds = (@set ds.G = I)
    dsθ = ds(θ)
    @unpack Cϕ,d = dsθ
    T = real(eltype(d))

    # compute approximate inverse ϕ Hessian used in gradient descent,
    # possibly from quadratic estimate
    if (Nϕ == :qe); Nϕ = quadratic_estimate(dsθ).Nϕ/2; end
    Hϕ⁻¹ = (Nϕ == nothing) ? Cϕ : pinv(pinv(Cϕ) + pinv(Nϕ))

    ϕ = (ϕstart != nothing) ? ϕstart : ϕ = zero(diag(Cϕ))
    tr = []
    state = nothing
    pbar = Progress(nsteps, (progress ? 0 : Inf), "MAP_marg: ")
    ProgressMeter.update!(pbar)

    for i=1:nsteps
        aggressive_gc && cuda_gc()
        g, state = δlnP_δϕ(
            ϕ, θ, ds;
            use_previous_MF = i>nsteps_with_meanfield_update,
            progress = false, return_state = true, previous_state = state,
            Nsims, Nbatch, weights, conjgrad_kwargs, aggressive_gc
        )
        ϕ += T(α) * Hϕ⁻¹ * g
        push!(tr, @dict(i,g,ϕ))
        next!(pbar, showvalues=[
            ("step",i),
            ("Ncg (data)", length(state.gQD.history)),
            ("Ncg (sims)", i<=nsteps_with_meanfield_update ? length(first(state.gQD_sims).history) : "0 (MF not updated)"),
            ("α",α)
        ])
    end

    set_distributed_dataset(nothing) # free memory, which got used inside δlnP_δϕ

    return ϕ, tr

end


text:
@doc raw"""
    stochastic_gradient_descent(M, gradF, x)

perform a stochastic gradient descent

# Input

* `M` a manifold ``\mathcal M``
* `gradF` – a gradient function, that either returns a vector of the subgradients
  or is a vector of gradients
* `x` – an initial value ``x ∈ \mathcal M``

# Optional
* `cost` – (`missing`) you can provide a cost function for example to track the function value
* `evaluation` – ([`AllocatingEvaluation`](@ref)) specify whether the gradient(s) works by
   allocation (default) form `gradF(M, x)` or [`MutatingEvaluation`](@ref) in place, i.e.
   is of the form `gradF!(M, X, x)` (elementwise).
* `evaluation_order` – (`:Random`) – whether
  to use a randomly permuted sequence (`:FixedRandom`), a per
  cycle permuted sequence (`:Linear`) or the default `:Random` one.
* `stopping_criterion` ([`StopAfterIteration`](@ref)`(1000)`)– a [`StoppingCriterion`](@ref)
* `stepsize` ([`ConstantStepsize`](@ref)`(1.0)`) a [`Stepsize`](@ref)
* `order_type` (`:RandomOder`) a type of ordering of gradient evaluations.
  values are `:RandomOrder`, a `:FixedPermutation`, `:LinearOrder`
* `order` - (`[1:n]`) the initial permutation, where `n` is the number of gradients in `gradF`.
* `retraction_method` – (`default_retraction_method(M)`) a `retraction(M,x,ξ)` to use.

# Output
* `x_opt` – the resulting (approximately critical) point of gradientDescent
OR
* `options` - the options returned by the solver (see `return_options`)
"""
function stochastic_gradient_descent(
    M::AbstractManifold, gradF::TDF, x; kwargs...
) where {TDF}
    x_res = allocate(x)
    copyto!(M, x_res, x)
    return stochastic_gradient_descent!(M, gradF, x_res; kwargs...)
end

@doc raw"""
    stochastic_gradient_descent!(M, gradF, x)

perform a stochastic gradient descent in place of `x`.

# Input

* `M` a manifold ``\mathcal M``
* `gradF` – a gradient function, that either returns a vector of the subgradients
  or is a vector of gradients
* `x` – an initial value ``x ∈ \mathcal M``

for all optional parameters, see [`stochastic_gradient_descent`](@ref).
"""
function stochastic_gradient_descent!(
    M::AbstractManifold,
    gradF::TDF,
    x;
    cost::TF=Missing(),
    direction::DirectionUpdateRule=StochasticGradient(zero_vector(M, x)),
    evaluation::AbstractEvaluationType=AllocatingEvaluation(),
    stopping_criterion::StoppingCriterion=StopAfterIteration(10000) |
                                          StopWhenGradientNormLess(1e-9),
    stepsize::Stepsize=ConstantStepsize(M),
    order_type::Symbol=:Random,
    order=collect(1:(gradF isa Function ? length(gradF(M, x)) : length(gradF))),
    retraction_method::AbstractRetractionMethod=default_retraction_method(M),
    return_options=false,
    kwargs...,
) where {TDF,TF}
    p = StochasticGradientProblem(M, gradF; cost=cost, evaluation=evaluation)
    o = StochasticGradientDescentOptions(
        x,
        zero_vector(M, x),
        direction;
        stopping_criterion=stopping_criterion,
        stepsize=stepsize,
        order_type=order_type,
        order=order,
        retraction_method=retraction_method,
    )
    o = decorate_options(o; kwargs...)
    resultO = solve(p, o)
    if return_options
        return resultO
    else
        return get_solver_result(resultO)
    end
end
function initialize_solver!(
    ::StochasticGradientProblem, o::StochasticGradientDescentOptions
)
    o.k = 1
    (o.order_type == :FixedRandom) && (shuffle!(o.order))
    return o
end
function step_solver!(
    p::StochasticGradientProblem, o::StochasticGradientDescentOptions, iter
)
    s, o.gradient = o.direction(p, o, iter)
    retract!(p.M, o.x, o.x, -s * o.gradient)
    o.k = ((o.k) % length(o.order)) + 1
    return o
end
get_solver_result(o::StochasticGradientDescentOptions) = o.x


text:
struct Simulation
  Φ::Problem
  nsims::Int #number of sims
  burnin::Int
  D::Vector{Float64} #dividend realizations
  V::Vector{Float64} #value realizations
  Y::Vector{Float64} #production
  Yₕ::Vector{Float64} #high state production
  Yₗ::Vector{Float64} #low state production
  K::Vector{Float64} #capital
end

function Simulation(Φ; nsims::Int = DEF_N_SIMS, burnin::Int)

  local s::Int = 1 #current state (index)
  local ls::Int #lagged state
  local iψ::Int = length(iψ) ÷ 2 #current policy state
  local liψ::Int #lagged policy state
  local ph::Float64

  local weights::Vector{Vector{Float64}} = (s::Int->vec(p[:,s])).(1:Φ.S)

  local V::Vector{Float64} = Vector{Float64}(undef, nsims) #value realizations
  local Y::Vector{Float64} = Vector{Float64}(undef, nsims) #current production
  local Yh::Vector{Float64} = Vector{Float64}(undef, nsims) #high state production
  local Yl::Vector{Float64} = Vector{Float64}(undef, nsims) #low state production
  local K::Vector{Float64} = Vector{Float64}(undef, nsims) #capital

  local samplefrom = (1,2,3)

  for i ∈ (-burnin):nsims
    #select the current state
    s = sample(samplefrom, weights[s])

    if i ≥ 1
      ph::Float64 =
      V[i] = Φ.F.V[iψ, s] #get the current value
      K[i] = ψ[1][iψ] #get the current capital level
      Yh[i] = ψ[2][iψ] #get the high state production
      Yl[i] = getyl(Yh[i],K[i], weights[s][1], 1 -weights[s][1], Φ.F) #derive low state production
      Y[i] = i==1 ? Yh[i] : Yl[i]
    end

    iψ = Φ.Π[iψ, s]
  end

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
struct AutoForwardDiff{chunksize} <: AbstractADType end
function AutoForwardDiff(chunksize=nothing)
    AutoForwardDiff{chunksize}()
end

struct AutoReverseDiff <: AbstractADType end
struct AutoTracker <: AbstractADType end
struct AutoZygote <: AbstractADType end

struct AutoFiniteDiff{T1,T2} <: AbstractADType
    fdtype::T1
    fdhtype::T2
end
AutoFiniteDiff(;fdtype = Val(:forward), fdhtype = Val(:hcentral)) =
                                                  AutoFiniteDiff(fdtype,fdhtype)

function default_chunk_size(len)
    if len < DEFAULT_CHUNK_THRESHOLD
        len
    else
        DEFAULT_CHUNK_THRESHOLD
    end
end

function instantiate_function(f, x, ::AbstractADType, p, num_cons = 0)
    grad   = f.grad   === nothing ? nothing : (G,x)->f.grad(G,x,p)
    hess   = f.hess   === nothing ? nothing : (H,x)->f.hess(H,x,p)
    hv     = f.hv     === nothing ? nothing : (H,x,v)->f.hv(H,x,v,p)
    cons   = f.cons   === nothing ? nothing : (x)->f.cons(x,p)
    cons_j = f.cons_j === nothing ? nothing : (res,x)->f.cons_j(res,x,p)
    cons_h = f.cons_h === nothing ? nothing : (res,x)->f.cons_h(res,x,p)

    OptimizationFunction{true,DiffEqBase.NoAD,typeof(f.f),typeof(grad),
                         typeof(hess),typeof(hv),typeof(cons),
                         typeof(cons_j),typeof(cons_h)}(f.f,
                         DiffEqBase.NoAD(),grad,hess,hv,cons,
                         cons_j,cons_h)
end

function instantiate_function(f::OptimizationFunction{true}, x, ::AutoForwardDiff{_chunksize}, p, num_cons = 0) where _chunksize

    chunksize = _chunksize === nothing ? default_chunk_size(length(x)) : _chunksize

    _f = (θ, args...) -> first(f.f(θ, p, args...))

    if f.grad === nothing
        gradcfg = (args...) -> ForwardDiff.GradientConfig(x -> _f(x, args...), x, ForwardDiff.Chunk{chunksize}())
        grad = (res, θ, args...) -> ForwardDiff.gradient!(res, x -> _f(x, args...), θ, gradcfg(args...), Val{false}())
    else
        grad = f.grad
    end

    if f.hess === nothing
        hesscfg = (args...) -> ForwardDiff.HessianConfig(x -> _f(x, args...), x, ForwardDiff.Chunk{chunksize}())
        hess = (res, θ, args...) -> ForwardDiff.hessian!(res, x -> _f(x, args...), θ, hesscfg(args...), Val{false}())
    else
        hess = f.hess
    end

    if f.hv === nothing
        hv = function (H,θ,v, args...)
            res = ArrayInterface.zeromatrix(θ)
            hess(res, θ, args...)
            H .= res*v
        end
    else
        hv = f.hv
    end

    if f.cons === nothing
        cons = nothing
        cons! = nothing
    else
        cons = θ -> f.cons(θ,p)
        cons! = (res, θ) -> (res .= f.cons(θ,p); res)
    end

    if cons !== nothing && f.cons_j === nothing
        cons_j = function (J, θ)
            cjconfig = ForwardDiff.JacobianConfig(cons, θ, ForwardDiff.Chunk{chunksize}())
            ForwardDiff.jacobian!(J, cons, θ, cjconfig)
        end
    else
        cons_j = f.cons_j
    end

    if cons !== nothing && f.cons_h === nothing
        cons_h = function (res, θ)
            for i in 1:num_cons
                hess_config_cache = ForwardDiff.HessianConfig(x -> cons(x)[i], θ,ForwardDiff.Chunk{chunksize}())
                ForwardDiff.hessian!(res[i], (x) -> cons(x)[i], θ, hess_config_cache,Val{false}())
            end
        end
    else
        cons_h = f.cons_h
    end

    return OptimizationFunction{true,AutoForwardDiff,typeof(f.f),typeof(grad),typeof(hess),typeof(hv),typeof(cons),typeof(cons_j),typeof(cons_h)}(f.f,AutoForwardDiff(),grad,hess,hv,cons,cons_j,cons_h)
end

function instantiate_function(f, x, ::AutoZygote, p, num_cons = 0)
    num_cons != 0 && error("AutoZygote does not currently support constraints")

    _f = (θ, args...) -> f(θ,p,args...)[1]
    if f.grad === nothing
        grad = (res, θ, args...) -> res isa DiffResults.DiffResult ? DiffResults.gradient!(res, Zygote.gradient(x -> _f(x, args...), θ)[1]) : res .= Zygote.gradient(x -> _f(x, args...), θ)[1]
    else
        grad = f.grad
    end

    if f.hess === nothing
        hess = function (res, θ, args...)
            if res isa DiffResults.DiffResult
                DiffResults.hessian!(res, ForwardDiff.jacobian(θ) do θ
                                                Zygote.gradient(x -> _f(x, args...), θ)[1]
                                            end)
            else
                res .=  ForwardDiff.jacobian(θ) do θ
                    Zygote.gradient(x ->_f(x, args...), θ)[1]
                end
            end
        end
    else
        hess = f.hess
    end

    if f.hv === nothing
        hv = function (H, θ, v, args...)
            _θ = ForwardDiff.Dual.(θ, v)
            res = DiffResults.GradientResult(_θ)
            grad(res, _θ, args...)
            H .= getindex.(ForwardDiff.partials.(DiffResults.gradient(res)),1)
        end
    else
        hv = f.hv
    end

    return OptimizationFunction{false,AutoZygote,typeof(f),typeof(grad),typeof(hess),typeof(hv),Nothing,Nothing,Nothing}(f,AutoZygote(),grad,hess,hv,nothing,nothing,nothing)
end

function instantiate_function(f, x, ::AutoReverseDiff, p=DiffEqBase.NullParameters(), num_cons = 0)
    num_cons != 0 && error("AutoReverseDiff does not currently support constraints")

    _f = (θ, args...) -> first(f.f(θ,p, args...))

    if f.grad === nothing
        grad = (res, θ, args...) -> ReverseDiff.gradient!(res, x -> _f(x, args...), θ, ReverseDiff.GradientConfig(θ))
    else
        grad = f.grad
    end

    if f.hess === nothing
        hess = function (res, θ, args...)
            if res isa DiffResults.DiffResult
                DiffResults.hessian!(res, ForwardDiff.jacobian(θ) do θ
                                                ReverseDiff.gradient(x -> _f(x, args...), θ)[1]
                                            end)
            else
                res .=  ForwardDiff.jacobian(θ) do θ
                    ReverseDiff.gradient(x ->_f(x, args...), θ)
                  end
            end
        end
    else
        hess = f.hess
    end

    if f.hv === nothing
        hv = function (H,θ,v, args...)
            _θ = ForwardDiff.Dual.(θ,v)
            res = DiffResults.GradientResult(_θ)
            grad(res, _θ, args...)
            H .= getindex.(ForwardDiff.partials.(DiffResults.gradient(res)),1)
        end
    else
        hv = f.hv
    end

    return OptimizationFunction{false,AutoReverseDiff,typeof(f),typeof(grad),typeof(hess),typeof(hv),Nothing,Nothing,Nothing}(f,AutoReverseDiff(),grad,hess,hv,nothing,nothing,nothing)
end

function instantiate_function(f, x, ::AutoTracker, p, num_cons = 0)
    num_cons != 0 && error("AutoTracker does not currently support constraints")
    _f = (θ, args...) -> first(f.f(θ, p, args...))

    if f.grad === nothing
        grad = (res, θ, args...) -> res isa DiffResults.DiffResult ? DiffResults.gradient!(res, Tracker.data(Tracker.gradient(x -> _f(x, args...), θ)[1])) : res .= Tracker.data(Tracker.gradient(x -> _f(x, args...), θ)[1])
    else
        grad = f.grad
    end

    if f.hess === nothing
        hess = (res, θ, args...) -> error("Hessian based methods not supported with Tracker backend, pass in the `hess` kwarg")
    else
        hess = f.hess
    end

    if f.hv === nothing
        hv = (res, θ, args...) -> error("Hessian based methods not supported with Tracker backend, pass in the `hess` and `hv` kwargs")
    else
        hv = f.hv
    end

    return OptimizationFunction{false,AutoTracker,typeof(f),typeof(grad),typeof(hess),typeof(hv),Nothing,Nothing,Nothing}(f,AutoTracker(),grad,hess,hv,nothing,nothing,nothing)
end

function instantiate_function(f, x, adtype::AutoFiniteDiff, p, num_cons = 0)
    num_cons != 0 && error("AutoFiniteDiff does not currently support constraints")
    _f = (θ, args...) -> first(f.f(θ, p, args...))

    if f.grad === nothing
        grad = (res, θ, args...) -> FiniteDiff.finite_difference_gradient!(res, x ->_f(x, args...), θ, FiniteDiff.GradientCache(res, x, adtype.fdtype))
    else
        grad = f.grad
    end

    if f.hess === nothing
        hess = (res, θ, args...) -> FiniteDiff.finite_difference_hessian!(res, x ->_f(x, args...), θ, FiniteDiff.HessianCache(x, adtype.fdhtype))
    else
        hess = f.hess
    end

    if f.hv === nothing
        hv = function (H, θ, v, args...)
            res = ArrayInterface.zeromatrix(θ)
            hess(res, θ, args...)
            H .= res*v
        end
    else
        hv = f.hv
    end

    return OptimizationFunction{false,AutoFiniteDiff,typeof(f),typeof(grad),typeof(hess),typeof(hv),Nothing,Nothing,Nothing}(f,adtype,grad,hess,hv,nothing,nothing,nothing)
end


text:
export SparseRowView, vals, cols

using SparseArrays
import SparseArrays: nnz
export nnz

"""
    SparseRowView(vals::AbstractArray{Data}, cols::AbstractArray{IndexType}, count::Integer=length(vals), start::Integer=1, stride::Integer=1)

Creates a view of a sparse row
"""
struct SparseRowView{Data, IndexType <: Integer}
    vals::AbstractArray{Data}
    cols::AbstractArray{IndexType}

    function SparseRowView(vals::AbstractArray{Data}, cols::AbstractArray{IndexType}
            ) where {Data, IndexType}
        if length(vals) != length(cols)
            throw(InvalidArgumentError("length(vals) = $(length(vals)) "
                    * "!= length(cols) = $(length(cols))"))
        end
        new{Data, IndexType}(vals, cols)
    end
end

function SparseRowView(vals::AbstractArray{Data}, cols::AbstractArray{IndexType},
        count::Integer, start::Integer=1, stride::Integer=1) where {Data, IndexType}
    SparseRowView(view(vals, range(start, step=stride, length=count)),
        view(cols, range(start, step=stride, length=count)))
end

function nnz(row::SparseRowView{Data, IndexType}) where{Data, IndexType}
    IndexType(length(row.vals))
end

"""
    vals(::SparseRowView{Data, IndexType})::AbstractArray{Data, 1}

Gets the values of the row
"""
vals(row::SparseRowView) = row.vals

"""
    cols(::SparseRowView{Data, IndexType})::AbstractArray{IndexType, 1}

Gets the column indices of the row
"""
cols(row::SparseRowView) = row.cols


text:
using Clang.Generators
using LMDB_jll

cd(@__DIR__)

include_dir = joinpath(LMDB_jll.artifact_dir, "include") |> normpath
lmdb_dir = include_dir

options = load_options(joinpath(@__DIR__, "generator.toml"))

# add compiler flags, e.g. "-DXXXXXXXXX"
args = get_default_args()
push!(args, "-I$include_dir")

headers = [joinpath(lmdb_dir, header) for header in readdir(lmdb_dir) if endswith(header, ".h")]
# there is also an experimental `detect_headers` function for auto-detecting top-level headers in the directory
# headers = detect_headers(clang_dir, args)

# create context
ctx = create_context(headers, args, options)

# run generator
build!(ctx, BUILDSTAGE_NO_PRINTING)

function rewrite!(e)
    if length(e)==1
        ex = first(e)
        if ex.head == :function && startswith(string(ex.args[1].args[1]),"mdb_")
            fname = ex.args[1].args[1]
            f2name = Symbol(string("_", fname))
            ex2 = Expr(:(=),deepcopy(ex.args[1]), Expr(:call, :checked_call, f2name,ex.args[1].args[2:end]...))
            ex.args[1].args[1] = f2name
            push!(e,ex2)
        end
    end
end

function rewrite!(dag::ExprDAG)
    for node in get_nodes(dag)
        exlist = get_exprs(node)
        rewrite!(exlist)
    end
end

rewrite!(ctx.dag)

# print
build!(ctx, BUILDSTAGE_PRINTING_ONLY)


text:
module Day18

using AdventOfCode2019

struct Key
    pos::CartesianIndex
    dist::Int
    doors::Array{Char,1}
    keys::Array{Char,1}
end

function day18(input::String = readInput(joinpath(@__DIR__, "..", "data", "day18.txt")))
    maze = parse_input(input)
    p1 = part1(maze)

    modify_maze!(maze)
    p2 = part2(maze)

    return [p1, p2]
end

function part1(maze::Array{Char,2})
    specialKeys = ['@']
    dists = generate_dists(maze, specialKeys)
    memo = Dict{Tuple{BitSet,BitSet,BitSet},Int}()
    keys = [x for x in maze if islowercase(x)]
    return shortest_distance!(dists, specialKeys, keys, Array{Char,1}(), memo)
end

function part2(maze::Array{Char,2})
    specialKeys = ['1', '2', '3', '4']
    dists = generate_dists(maze, specialKeys)
    memo = Dict{Tuple{BitSet,BitSet,BitSet},Int}()
    keys = [x for x in maze if islowercase(x)]
    return shortest_distance!(dists, specialKeys, keys, Array{Char,1}(), memo)
end

function parse_input(input::String)
    m = hcat(split.(split(strip(input), '\n'), "")...)
    maze = Array{Char,2}(undef, size(m, 2), size(m, 1))
    for i = 1:size(m, 2)
        for j = 1:size(m, 1)
            maze[i, j] = Char(m[j, i][1])
        end
    end
    return maze
end

function generate_dists(maze, specialKeys)
    dists = Dict{Char,Dict{Char,Key}}()
    keys = [x for x in maze if islowercase(x)]
    for k in union(keys, specialKeys)
        dists[k] = flood(maze, k)
    end
    return dists
end

function flood(maze::Array{Char,2}, start::Char)
    cmaze = copy(maze)
    startInd = findfirst(x -> x == start, maze)
    distDict = Dict{Char,Key}()
    queue = Array{Key,1}()
    push!(queue, Key(startInd, 0, [], []))
    while length(queue) > 0
        currentKey = popfirst!(queue)
        if islowercase(cmaze[currentKey.pos]) || cmaze[currentKey.pos] in ('@', '1', '2', '3', '4')
            distDict[cmaze[currentKey.pos]] = deepcopy(currentKey)
            if islowercase(cmaze[currentKey.pos])
                push!(currentKey.keys, cmaze[currentKey.pos])
            end
        elseif isuppercase(cmaze[currentKey.pos])
            push!(currentKey.doors, cmaze[currentKey.pos])
        end
        cmaze[currentKey.pos] = '#'
        for (k, l) in ((-1, 0), (1, 0), (0, -1), (0, 1))
            if cmaze[(currentKey.pos.I .+ (k, l))...] != '#'
                push!(
                    queue, Key(CartesianIndex(currentKey.pos.I .+ (k, l)),
                    currentKey.dist + 1, copy(currentKey.doors), copy(currentKey.keys)))
            end
        end
    end
    return distDict
end

function modify_maze!(maze::Array{Char,2})
    atPos = findfirst(x -> x == '@', maze)
    maze[atPos[1]-1:atPos[1]+1, atPos[2]-1:atPos[2]+1] .= ['1' '#' '2';
                                                           '#' '#' '#';
                                                           '3' '#' '4']
end

function shortest_distance!(dists::Dict{Char,Dict{Char,Key}}, current::Array{Char}, remainingKeys::Array{Char,1}, unlockedDoors::Array{Char,1}, memo::Dict{Tuple{BitSet,BitSet,BitSet},Int})
    length(remainingKeys) == 0 && return 0
    candids = [Array{Char,1}() for i = 1:length(current)]
    for i = 1:length(current)
        for key in remainingKeys
            haskey(dists[current[i]], key) || continue
            iscandid = true
            # Robot must be able to reach `key`. Otherwise it cannot
            # be considered as a candidate.
            for door in dists[current[i]][key].doors
                if !(door in unlockedDoors)
                    iscandid = false
                    break
                end
            end
            !iscandid && continue
            # If the robot would collect other keys on its way from
            # `current` to `key`, we do not consider `key` to be a
            # candidate.
            for k in dists[current[i]][key].keys
                if k in remainingKeys
                    iscandid = false
                    break
                end
            end
            if iscandid
                push!(candids[i], key)
            end
        end
    end
    state = (BitSet(Int.(current)), BitSet(Int.(remainingKeys)), BitSet(Int.(unlockedDoors)))
    haskey(memo, state) && return memo[state]
    total = minimum(
        [dists[current[i]][c].dist + shortest_distance!(
            dists, replace(x -> x == current[i] ? c : x, current), setdiff(remainingKeys, c), union(unlockedDoors, uppercase.(dists[current[i]][c].keys), [uppercase(c)]), memo) for i=1:length(current) for c in candids[i]])
    memo[state] = total
    return total
end

end # module


text:
using {{{PKG}}}
using Test

@testset "{{{PKG}}}.jl" begin
    # Write your tests here.
end


text:
function can_move(cur_player::AbstractPlayer)
  cur_player.game.in_focus || return false
  cur_player.state.is_playing || return false

  return true
end


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

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 el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ eceac3a0-1e4a-11eb-19d7-e1a6af47842e
begin
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.add(["DifferentialEquations", "Plots", "PlutoUI", "Distributions"])
	using DifferentialEquations
	using PlutoUI
	using Plots
	using Distributions
end

# ╔═╡ f880960c-1eb5-11eb-08d4-656728e0297e
md"_*reading notes from 《METABOLIC GROWTH THEORY》, book author: Chen ping.*_"

# ╔═╡ 506c8354-1e42-11eb-3537-d9c055da2552
md"## Problems?"

# ╔═╡ e1ee8568-1e40-11eb-00d2-3d279f8e9bb3
md"_**What is science?**_"

# ╔═╡ d1924e1a-1e41-11eb-3640-ab6e3784c8ec
md"> Science is an institution, it is a method, it is an accumulated tradition of kownledge, it is a major factor in the development of production, and it is a powerful force in the information of the universe"

# ╔═╡ 59efff32-1e42-11eb-224f-bd5a0929bd69
md"_**China led the world in science, but lagged behind only after the Ming and Qing Dynasties?**_"

# ╔═╡ a8ca358c-1e42-11eb-0122-4dedb09ca5e1
md"> so we say behind only after the Ming and Qing Dynasties, we are behind the Western. it will reduce the stigma, but it impedes our understanding of the weaknesses of traditional chinese secience.
> China was always backward in the scientific systerm and the foundation of scientific method, and the independent and complete scientific system never emerged from the feudal society.Not to admit the backwardness of system and method was just the staring point of the backwardness of modern china."

# ╔═╡ e13f067a-1e44-11eb-044f-e5dc23a70983
md"_**What is Technology?**_"

# ╔═╡ 04f3d500-1e45-11eb-2ed6-1b1057a090b8
md"> Technology is the use of scientific methods and knowledge to solve the problems facing human life in order to obtain the greatest possible economic benefits."

# ╔═╡ 704a389a-1e46-11eb-2717-f17fc35381e0
md"# Content"

# ╔═╡ 867332b4-1e46-11eb-1e49-470a9cdded52
md"## The Origin of division of Labor and a stochastic model of social differentiation"

# ╔═╡ 456b046c-1e47-11eb-0aea-510bda71e16e
md"""
We note that **_cultural factors_** played an important role in the origins of capitalism and modern science.

he proposed a one-dimensional model of individualism, assuming a horizontal axis, with highly individualistic Europe and the United States at one extreme and societies with rigid division of Labour like bees at the other.(Kikuchi, 1981)

**_Information diffusion model_**:Use Logistic equation to sociology.

$\frac{\partial n}{\partial t} = kn(N-n) - dn(1 - \alpha \frac{n}{N})$

noted:$N$ represents the population size, $n$ is the number of people who have mastered the information, $(N-n)$ is the number of learners who need to receive new information, and $k$ is the growth rate of the number of people who receive new information in the learning process. The Last term represents the elimination rate or forgetting rate. $\alpha$ is a measure of how repulsive people are to unfamiliar things
"""

# ╔═╡ 522f7894-1e73-11eb-0131-65336d6da211
md"""

so $N$ : $(@bind N Slider(10:20; default=10, show_value=true))

const $k$ : $(@bind k Slider(0.00 : 0.01 : 0.50; default=0.02, show_value=true))

rate α : $(@bind α Slider(-0.50 : 0.01 : 0.50; default=0.00, show_value=true))
"""

# ╔═╡ 39d1607e-1e74-11eb-2c9c-ab49addfdc1d
begin
	Infordiff(u,p,t)=k*u*(N-u) - 0.02 * u*(1.0 - α* u/N)
	u0 = 0.1
	tspan = (0.0,100.0)
	prob = ODEProblem(Infordiff,u0,tspan)
	sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
	plot(sol)
end

# ╔═╡ b6bce1e4-1e74-11eb-2ba5-8f641457daa8
function Infor(ina)
	plot(location=4,legend = :outertopright)
	for i in ina
		Infordiff1(u,p,t)=0.02*u*(10-u) - 0.1u * (1 - i * u/10)
		u0 = 0.1
		tspan = (0.0,100.0)
		prob = ODEProblem(Infordiff1,u0,tspan)
		sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
		plot!(sol, label="α = $(i)")
	end
	plot!(axis=false)
end

# ╔═╡ 421f5f96-1e75-11eb-0783-251da9b8ca47
Infor(-0.9:0.2:0.9)

# ╔═╡ 0f29b070-1eb4-11eb-28bb-d723b952495e
md"""we can solve that if the $-1<\alpha <0$, $\alpha$ is a measure of how repulsive people are to unfamiliar things. use it to measure society's preference for risky behavior, more closely to $-1$, you are more like risky or you can more receptive to new thing.

if $\alpha$ in $[0,1]$,we measure pepole's rejection of unfamiliar things.

from the plot above, we can see the number of equilibrium from $-0.9$ to $0.9$,the $n$ is different,the smaller $\alpha$ will get smaller equilibrium $n$.

**_Conclusion_: $n_{\alpha <0} < n_{\alpha=0} < n_{\alpha>0}$, The balenced population or market size of the individualist culture is smaller than that of the collectivist culture.**
"""

# ╔═╡ fb3c4b76-1eb4-11eb-2824-4104209213b2
md" ## A fluctuating environment"

# ╔═╡ a2ac3c50-1eb6-11eb-17cd-359e9d98b756
md"if you are dealing with a fluctuating environment, use this function, here:"

# ╔═╡ Cell order:
# ╟─f880960c-1eb5-11eb-08d4-656728e0297e
# ╟─506c8354-1e42-11eb-3537-d9c055da2552
# ╟─e1ee8568-1e40-11eb-00d2-3d279f8e9bb3
# ╟─d1924e1a-1e41-11eb-3640-ab6e3784c8ec
# ╟─59efff32-1e42-11eb-224f-bd5a0929bd69
# ╟─a8ca358c-1e42-11eb-0122-4dedb09ca5e1
# ╟─e13f067a-1e44-11eb-044f-e5dc23a70983
# ╟─04f3d500-1e45-11eb-2ed6-1b1057a090b8
# ╟─704a389a-1e46-11eb-2717-f17fc35381e0
# ╟─867332b4-1e46-11eb-1e49-470a9cdded52
# ╟─456b046c-1e47-11eb-0aea-510bda71e16e
# ╟─eceac3a0-1e4a-11eb-19d7-e1a6af47842e
# ╟─522f7894-1e73-11eb-0131-65336d6da211
# ╟─39d1607e-1e74-11eb-2c9c-ab49addfdc1d
# ╟─b6bce1e4-1e74-11eb-2ba5-8f641457daa8
# ╟─421f5f96-1e75-11eb-0783-251da9b8ca47
# ╟─0f29b070-1eb4-11eb-28bb-d723b952495e
# ╟─fb3c4b76-1eb4-11eb-2824-4104209213b2
# ╟─a2ac3c50-1eb6-11eb-17cd-359e9d98b756

max_stars_repo_path: Visualization/information diff.jl

max_stars_repo_name: ZhouZhuofei/juliaDL

max_issues_repo_path: Visualization/information diff.jl

max_issues_repo_name: ZhouZhuofei/juliaDL

max_forks_repo_path: Visualization/information diff.jl

max_forks_repo_name: ZhouZhuofei/juliaDL


text:
module MNIST_BATCH
using Flux
using Flux.Data.MNIST, Statistics
using Flux: onehotbatch, onecold, crossentropy,throttle
using Base.Iterators: repeated,partition

using CUDAnative
using CuArrays
CuArrays.allowscalar(false)

#=
Very important !!
ϵ is used to prevent loss NaN
=#
const ϵ = 1.0f-32

# Load training labels and images from Flux.Data.MNIST
@info("Loading data...")

train_imgs = MNIST.images()
train_labels = MNIST.labels()

# use 1nd GPU : default
CUDAnative.device!(0)
# use 2nd GPU
#CUDAnative.device!(1)

# Bundle images together with labels and group into minibatch
function make_minibatch(imgs,labels,batch_size)
  X = hcat(float.(reshape.(imgs,:))...) |> gpu
  Y = float.(onehotbatch(labels,0:9)) |> gpu

  data_set = [(X[:,i],Y[:,i]) for i in partition(1:length(labels),batch_size)]
  return data_set
end

@info("Making model...")
# Model
m = Chain(
  Dense(28^2,32,relu),
  Dense(32,10),
  softmax
) |> gpu
loss(x,y) = crossentropy(m(x) .+ ϵ, y)
accuracy(x,y) = mean(onecold(m(x)|>cpu) .== onecold(y|>cpu))

batch_size = 500
train_dataset = make_minibatch(train_imgs,train_labels,batch_size)

opt = ADAM()

@info("Training model...")

epochs = 200
# used for plots
accs = Array{Float32}(undef,0)

dataset_len = length(train_dataset)
for i in 1:epochs
  for (idx,dataset) in enumerate(train_dataset)
    Flux.train!(loss,params(m),[dataset],opt)
    acc = accuracy(dataset...)
    if idx == dataset_len
      @info("Epoch# $(i)/$(epochs) - loss: $(loss(dataset...)), accuracy: $(acc)")
      push!(accs,acc)
    end
  end
end

# Test Accuracy
tX = hcat(float.(reshape.(MNIST.images(:test),:))...) |> gpu
tY = float.(onehotbatch(MNIST.labels(:test),0:9)) |> gpu

println("Test loss:", loss(tX,tY))
println("Test accuracy:", accuracy(tX,tY))

end

using Plots;gr()
plot(MNIST_BATCH.accs)

max_stars_repo_name: mrchaos/model-zoo

max_issues_repo_name: mrchaos/model-zoo

max_forks_repo_name: mrchaos/model-zoo


text:
# Implements the selection sort algorithm
using Debug

function swap2!(s::Array, i::Integer, k::Integer)
  temp = s[i]
  s[i] = s[k]
  s[k] = temp
end

function omerselect(s::Array)
	for i=1:length(s)
		minIndex = i
		for k=i+1:length(s)
	  		(s[minIndex] > s[k]) && (minIndex = k)
		end
		(i!=minIndex) && (swap2!(s, minIndex, i))
	end
	return s
end

function test()
  a = [8, 4, 2, 5, 42]
  b = rand(4,1)
  println(omerselect(a))

  println("=======")
  omerselect(b)
end
test()

max_stars_repo_name: kirnap/Machine-Learning-in-Julia

max_issues_repo_name: kirnap/Machine-Learning-in-Julia

max_forks_repo_name: kirnap/Machine-Learning-in-Julia


text: function f () 1 end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
module Graphs
import Base.add

abstract Value

type Node
    value::Value
    refs::Set{Node}
    replacement::Union{Node,Nothing}

    Node(value) = new(value, Set{Node}(), nothing)
end

type Graph
    nodes::Dict{Value, Node}

    Graph() = new(Dict{Value,Node}())
end

add(g::Graph, v::Value) = (has(g.nodes, v) ? g.nodes[v] : g.nodes[v] = Node(v))

function subs!(g::Graph, from::Value, to::Value)
    has(g.nodes, from) ? subs!(g, g.nodes[from], to) : nothing
end

function subs!(g::Graph, from::Node, to_value::Value)
    if has(g.nodes, to_value)
        subs!(g, from, g.nodes[to_value])
    else
        from.value = to_value
        g.nodes[to_value] = from
        return
    end
end

function subs!(g::Graph, from::Node, to::Node)
    from.replacement = to
    g.nodes[from.value] = to
    # replace from with to in all references to from
    for refnode in from.refs
        subs!(g, refnode, subs(refnode.value, from, to))
    end
    add_each(to.refs, from.refs)
end

end # module


text:
function validtreeconstraints(
    pd::PopulationData,
    bt::BinaryTree,
    tree::JuMP.Model,
    assign::JuMP.JuMPArray{JuMP.Variable},
    outgroupnode::Int,
    nlevels::Int,
    nmixtures::Int)

    npop = pd.npop
    levels = 1:nlevels
    leaves = getleaves(bt)

    # each population gets assigned levels
    JuMP.@constraint(tree,
        [a=1:npop,l=levels],
        sum(assign[a,leaves,l]) == 1)
    # each node gets at most one population
    if nlevels > 1
        JuMP.@variable(tree, assign1[1:npop,leaves] >= 0)
        JuMP.@constraint(tree,
            [a=1:npop,n=leaves,l=levels],
            assign1[a,n] >= assign[a,n,l])
        JuMP.@constraint(tree,
            [n=leaves],
            sum(assign1[a,n] for a in 1:npop) <= 1)
        JuMP.@constraint(tree,
            sum(assign1) <= nmixtures)
    else
        JuMP.@constraint(tree,
            [n=leaves],
            sum(assign[1:npop,n,1]) <= 1)
    end
    # assign outgroup to a node
    JuMP.@constraint(tree,
        [l=levels],
        assign[pd.outgroup,outgroupnode,l] == 1)

end

# Juan Pablo trickery to improve branch-and-bound performance
function binaryencodingconstraints(
    pd::PopulationData,
    bt::BinaryTree,
    tree::JuMP.Model,
    assign::JuMP.JuMPArray{JuMP.Variable})

    JuMP.@variable(tree, codeselect[1:pd.npop,1:bt.depth], Bin)
    JuMP.@constraint(tree,
        [a=1:pd.npop,m=1:bt.depth],
        codeselect[a,m] == sum(bt.codes[u][m]*assign[a,u,1] for u in getleaves(bt)))

end

function warmstartunmixed(tp::Union{NodeTreeProblem,TreeProblem};
    timelimit::Int = 30,
    solver = Gurobi.GurobiSolver(OutputFlag = 0, TimeLimit = timelimit))
    @assert tp.nlevels > 1
    tp0 = TreeProblem(tp.pd, tp.bt, solver = solver)
    JuMP.solve(tp0.model)
    solution0 = JuMP.getvalue(tp0.assign);
    for a in 1:tp.pd.npop, u in getnodes(tp.bt, tp.bt.depth), l in 1:tp.nlevels
        JuMP.setvalue(tp.assign[a,u,l], JuMP.getvalue(tp0.assign[a,u,1]))
        JuMP.setvalue(tp.model[:assign1][a,u], JuMP.getvalue(tp0.assign[a,u,1]))
    end
end

# fixes population a to be un-admixed
function unmix(tp::Union{NodeTreeProblem,TreeProblem}, a::Int)
    JuMP.@constraint(tp.model,
        [n=getleaves(tp.bt),l=2:tp.nlevels],
        tp.assign[a,n,l] == tp.assign[a,n,l-1])
end

# fixes pop to be un-admixed
function unmix(tp::Union{NodeTreeProblem,TreeProblem}, pop::String)
    a = findfirst(tp.pd.pops .== pop)
    if a == nothing
        @error("Invalid population $(pop)")
    else
        unmix(tp, a)
    end
end

# fixes population a to node u
function fix(tp::Union{NodeTreeProblem,TreeProblem}, a::Int, u::Int)
    JuMP.@constraint(tp.model,
        [l=1:tp.nlevels],
        tp.assign[a,u,l] == 1)
end

# fixes pop to node u
function fix(tp::Union{NodeTreeProblem,TreeProblem}, pop::String, u::Int)
    a = findfirst(tp.pd.pops .== pop)
    if a == nothing
        @error("Invalid population $(pop)")
    else
        fix(tp, a, u)
    end
end

function removesolution(
    tp::Union{NodeTreeProblem,TreeProblem},
    solution::JuMP.JuMPArray{Float64})
    expr = 0
    for a in 1:tp.pd.npop, n in getnodes(tp.bt, tp.bt.depth), l in 1:tp.nlevels
        if solution[a,n,l] > 0
            expr += tp.assign[a,n,l]
        else
            expr += (1-tp.assign[a,n,l])
        end
    end
    JuMP.@constraint(tp.model, expr <= JuMP.getvalue(expr)-1)
end