text:
module ActifMPCCmod

import Relaxation
import ParamSetmod
import NLPModels

"""
Type MPCC_actif : problème MPCC pénalisé avec slack et ensemble des contraintes actives

liste des constructeurs :
MPCC_actif(nlp::NLPModels.AbstractNLPModel,r::Float64,s::Float64,t::Float64,nb_comp::Int64,paramset::ParamSetmod.ParamSet,direction::Function,linesearch::Function)
MPCC_actif(nlp::NLPModels.AbstractNLPModel,r::Float64,s::Float64,t::Float64,w::Any,paramset::ParamSetmod.ParamSet,direction::Function,linesearch::Function)

liste des méthodes :
updatew(ma::MPCC_actif)
setf(ma::MPCC_actif, f::Function, xj::Any)
setw(ma::MPCC_actif, w::Any)
setbeta(ma::MPCC_actif,b::Float64)
sethess(ma::MPCC_actif,Hess::Array{Float64,2})

liste des accesseurs :
evalx(ma::MPCC_actif,x::Vector)
evald(ma::MPCC_actif,d::Vector)
redd(ma::MPCC_actif,d::Vector)
obj(ma::MPCC_actif,x::Vector)
ExtddDirection(ma::ActifMPCCmod.MPCC_actif,dr::Vector,xp::Vector,step::Float64)
grad(ma::MPCC_actif,x::Vector)
grad(ma::MPCC_actif,x::Vector,gradf::Vector)
hess(ma::MPCC_actif,x::Vector)
hess(ma::MPCC_actif,x::Vector,H::Array{Float64,2})

liste des fonctions :
LSQComputationMultiplier(ma::MPCC_actif,gradpen::Vector,xj::Vector)
RelaxationRule(ma::ActifMPCCmod.MPCC_actif,xj::Vector,lg::Vector,lh::Vector,lphi::Vector,wmax::Any)
PasMaxComp(ma::MPCC_actif,x::Vector,d::Vector)
AlphaChoix(alpha::Float64,alpha1::Float64,alpha2::Float64)
AlphaChoix(alpha::Float64,alpha1::Float64,alpha2::Float64,alpha3::Float64,alpha4::Float64)
PasMaxBound(ma::MPCC_actif,x::Vector,d::Vector) #TO DO
PasMax(ma::MPCC_actif,x::Vector,d::Vector)
"""

#TO DO List :
#Major :
#- Changer les w en set ou en liste
#- intégrer les contraintes de bornes sur x dans PasMax
#Minor :
#- copie de code dans PasMax
#- Plutôt que de passer le point courant -> stocker en point initial de mpcc ?

type MPCC_actif
 nlp::NLPModels.AbstractNLPModel # en fait on demande juste une fonction objectif, point initial, contraintes de bornes
 r::Float64
 s::Float64
 t::Float64

 w::Array{Bool,2} #matrice à 2 colonnes et de longueur 2nb_comp -- sparse

 n::Int64 #dans le fond est optionnel si on a nb_comp
 nb_comp::Int64

 #ensembles d'indices :
 w1::Array{Int64,1} # ensemble des indices (entre 0 et nb_comp) où la contrainte yG>=-r est active
 w2::Array{Int64,1} # ensemble des indices (entre 0 et nb_comp) où la contrainte yH>=-r est active
 w3::Array{Int64,1} # ensemble des indices (entre 0 et nb_comp) où la contrainte yG<=s+t*theta(yH,r) est active
 w4::Array{Int64,1} # ensemble des indices (entre 0 et nb_comp) où la contrainte yH<=s+t*theta(yG,r) est active
 wcomp::Array{Int64,1} #ensemble des indices (entre 0 et nb_comp) où la contrainte Phi<=0 est active
 w13c::Array{Int64,1} #ensemble des indices où les variables yG sont libres
 w24c::Array{Int64,1} #ensemble des indices où les variables yH sont libres
 wc::Array{Int64,1} #ensemble des indices des contraintes où yG et yH sont libres
 wcc::Array{Int64,1} #ensemble des indices des contraintes où yG et yH sont fixés

 wnew::Array{Bool,2} #dernières contraintes ajoutés

 #paramètres pour le calcul de la direction de descente
 crho::Float64 #constant such that : ||c(x)||_2 \approx crho*rho
 beta::Float64 #paramètre pour gradient conjugué
 Hess::Array{Float64,2} #inverse matrice hessienne approximée
 #Hd::Vector #produit inverse matrice hessienne et gradient (au lieu de la hessienne entière)

 paramset::ParamSetmod.ParamSet
 direction::Function #fonction qui calcul la direction de descente
 linesearch::Function #fonction qui calcul la recherche linéaire
end

"""
Constructeur recommandé pour MPCC_actif
"""
function MPCC_actif(nlp::NLPModels.AbstractNLPModel,r::Float64,s::Float64,t::Float64,nb_comp::Int64,paramset::ParamSetmod.ParamSet,direction::Function,linesearch::Function)

  n=length(nlp.meta.x0)-2*nb_comp
  xk=nlp.meta.x0[1:n]
  ygk=nlp.meta.x0[n+1:n+nb_comp]
  yhk=nlp.meta.x0[n+nb_comp+1:n+2*nb_comp]

  w=zeros(Bool,2*nb_comp,2)

  for l=1:nb_comp
   if ygk[l]==nlp.meta.lvar[n+l]
    w[l,1]=true;
   elseif ygk[l]==Relaxation.psi(yhk[l],r,s,t)
    w[l+nb_comp,1]=true;
   end
   if yhk[l]==nlp.meta.lvar[n+l+nb_comp]
    w[l,2]=true;
   elseif yhk[l]==Relaxation.psi(ygk[l],r,s,t)
    w[l+nb_comp,2]=true;
   end
  end

 return MPCC_actif(nlp,r,s,t,w,paramset,direction,linesearch)
end

function MPCC_actif(nlp::NLPModels.AbstractNLPModel,r::Float64,s::Float64,t::Float64,w::Array{Bool,2},paramset::ParamSetmod.ParamSet,direction::Function,linesearch::Function)

 nb_comp=Int(size(w,1)/2)
 n=length(nlp.meta.x0)-2*nb_comp
 w1=find(w[1:nb_comp,1])
 w2=find(w[1:nb_comp,2])
 w3=find(w[nb_comp+1:2*nb_comp,1])
 w4=find(w[nb_comp+1:2*nb_comp,2])
 wcomp=find(w[nb_comp+1:2*nb_comp,1] .| w[nb_comp+1:2*nb_comp,2])
 w13c=find(.!w[1:nb_comp,1] .& .!w[nb_comp+1:2*nb_comp,1])
 w24c=find(.!w[1:nb_comp,2] .& .!w[nb_comp+1:2*nb_comp,2])
 wc=find(.!w[1:nb_comp,1] .& .!w[1:nb_comp,2] .& .!w[nb_comp+1:2*nb_comp,1] .& .!w[nb_comp+1:2*nb_comp,2])
 wcc=find((w[1:nb_comp,1] .| w[nb_comp+1:2*nb_comp,1]) .& (w[1:nb_comp,2] .| w[nb_comp+1:2*nb_comp,2]))

 wnew=zeros(Bool,0,0)

 crho=1.0
 beta=0.0
 Hess=eye(n+2*nb_comp)

 return MPCC_actif(nlp,r,s,t,w,n,nb_comp,w1,w2,w3,w4,wcomp,w13c,w24c,wc,wcc,wnew,crho,beta,Hess,paramset,direction,linesearch)
end

"""
Methodes pour le MPCC_actif
"""
#Mise à jour des composantes liés à w
function updatew(ma::MPCC_actif)

 ma.w1=find(ma.w[1:ma.nb_comp,1])
 ma.w2=find(ma.w[1:ma.nb_comp,2])
 ma.w3=find(ma.w[ma.nb_comp+1:2*ma.nb_comp,1])
 ma.w4=find(ma.w[ma.nb_comp+1:2*ma.nb_comp,2])
 ma.wcomp=find(ma.w[ma.nb_comp+1:2*ma.nb_comp,1] .| ma.w[ma.nb_comp+1:2*ma.nb_comp,2])
 ma.w13c=find(.!ma.w[1:ma.nb_comp,1] .& .!ma.w[ma.nb_comp+1:2*ma.nb_comp,1])
 ma.w24c=find(.!ma.w[1:ma.nb_comp,2] .& .!ma.w[ma.nb_comp+1:2*ma.nb_comp,2])
 ma.wc=find(.!ma.w[1:ma.nb_comp,1] .& .!ma.w[1:ma.nb_comp,2] .& .!ma.w[ma.nb_comp+1:2*ma.nb_comp,1] .& .!ma.w[ma.nb_comp+1:2*ma.nb_comp,2])
 ma.wcc=find((ma.w[1:ma.nb_comp,1] .| ma.w[ma.nb_comp+1:2*ma.nb_comp,1]) .& (ma.w[1:ma.nb_comp,2] .| ma.w[ma.nb_comp+1:2*ma.nb_comp,2]))
 return ma
end

#Mise à jour de w
function setw(ma::MPCC_actif, w::Array{Bool,2})

 ma.wnew=w .& .!ma.w
 ma.w=w
 return updatew(ma)
end

function setbeta(ma::MPCC_actif,b::Float64)
 ma.beta=b
 return ma
end

function setcrho(ma::MPCC_actif,crho::Float64)
 ma.crho=crho
 return ma
end

function sethess(ma::MPCC_actif,Hess::Array{Float64,2})
 ma.Hess=Hess
 return ma
end

"""
Renvoie le vecteur x=[x,yg,yh] au complet
"""
function evalx(ma::MPCC_actif,x::Vector)
 #construction du vecteur de taille n+2nb_comp que l'on évalue :
 xf=ma.s*ones(ma.n+2*ma.nb_comp)
 xf[1:ma.n]=x[1:ma.n]
 xf[ma.w13c+ma.n]=x[ma.n+1:ma.n+length(ma.w13c)]
 xf[ma.w24c+ma.n+ma.nb_comp]=x[ma.n+length(ma.w13c)+1:ma.n+length(ma.w13c)+length(ma.w24c)]

 #on regarde les variables yG fixées :
 xf[ma.w1+ma.n]=ma.nlp.meta.lvar[ma.w1+ma.n]
 xf[ma.w3+ma.n]=Relaxation.psi(xf[ma.w3+ma.n+ma.nb_comp],ma.r,ma.s,ma.t)
 #on regarde les variables yH fixées :
 xf[ma.w2+ma.n+ma.nb_comp]=ma.nlp.meta.lvar[ma.w2+ma.n+ma.nb_comp]
 xf[ma.w4+ma.n+ma.nb_comp]=Relaxation.psi(xf[ma.w4+ma.n],ma.r,ma.s,ma.t)

 return xf
end

"""
Renvoie la direction d au complet (avec des 0 aux actifs)
"""
function evald(ma::MPCC_actif,d::Vector)
 df=zeros(ma.n+2*ma.nb_comp)
 df[1:ma.n]=d[1:ma.n]
 df[ma.w13c+ma.n]=d[ma.n+1:ma.n+length(ma.w13c)]
 df[ma.w24c+ma.nb_comp+ma.n]=d[ma.n+length(ma.w13c)+1:ma.n+length(ma.w13c)+length(ma.w24c)]
 return df
end

"""
Renvoie la direction d réduite
"""
function redd(ma::MPCC_actif,d::Vector)

  df=zeros(ma.n+length(ma.w13c)+length(ma.w24c))
  df[1:ma.n]=d[1:ma.n]
  df[ma.n+1:ma.n+length(ma.w13c)]=d[ma.w13c+ma.n]
  df[ma.n+length(ma.w13c)+1:ma.n+length(ma.w13c)+length(ma.w24c)]=d[ma.w24c+ma.nb_comp+ma.n]

 return df
end

function redd(ma::MPCC_actif,d::Vector,w::Array{Int64,1})

  df=zeros(length(w))
  df[1:length(w)]=d[w]

 return df
end

"""
Evalue la fonction objectif d'un MPCC actif : x
"""
function obj(ma::MPCC_actif,x::Vector)

 if length(x)==ma.n+2*ma.nb_comp
  return NLPModels.obj(ma.nlp,x)
 else
  return NLPModels.obj(ma.nlp,evalx(ma,x))
 end

end

"""
Calcul la direction étendue à partir de la direction du sous-domaine
complète une version étendue de la direction de descente :
dr : la direction dans le sous-espace actif
xp : le nouvel itéré
step : le pas utilisé pour calculé xp

utilise la formule suivante :
yG+=yG+alpha*dyG --> dyG=(yG+-yG)/alpha
"""
function ExtddDirection(ma::ActifMPCCmod.MPCC_actif,dr::Vector,xp::Vector,step::Float64)
 d=evald(ma,dr) #evald rempli les trous par des 0
 x=evalx(ma,xp)

 d[1:ma.n]=dr[1:ma.n]

 psip=Relaxation.psi(x[ma.n+1:ma.n+2*ma.nb_comp],ma.r,ma.s,ma.t)
 psi=Relaxation.psi(x[ma.n+1:ma.n+2*ma.nb_comp]-step*d[ma.n+1:ma.n+2*ma.nb_comp],ma.r,ma.s,ma.t)
 #d[ma.w1]=0 #yG fixé
 #d[ma.w2]=0 #yH fixé

 d[ma.n+ma.w3]=(psip[ma.w3+ma.nb_comp]-psi[ma.w3+ma.nb_comp])/step #yG fixé
 d[ma.n+ma.nb_comp+ma.w4]=(psip[ma.w4]-psi[ma.w4])/step #yH fixé
 #d[ma.n+ma.w13c]=dr[ma.n+ma.w13c] #yG est libre
 #d[ma.n+ma.nb_comp+ma.w24c]=dr[ma.n+ma.nb_comp+ma.w24c] #yG est libre

 return d
end

"""
Evalue le gradient de la fonction objectif d'un MPCC actif
x est le vecteur réduit
"""
function grad(ma::MPCC_actif,x::Vector)

 #on calcul xf le vecteur complet
 xf=evalx(ma,x)
 #construction du vecteur gradient de taille n+2nb_comp
 gradf=NLPModels.grad(ma.nlp,xf)

 return length(x)==ma.n+2*ma.nb_comp?gradf:grad(ma,x,gradf)
end

function grad(ma::MPCC_actif,x::Vector,gradf::Vector)

 #on calcul xf le vecteur complet
 xf=evalx(ma,x)
 #construction du vecteur gradient de taille n+2nb_comp
 #gradf=NLPModels.grad(ma.nlp,xf)

 gradg=Array{Float64}
 # Conditionnelles pour gérer le cas où w1 et w3 est vide
 #if isempty(ma.w1) && isempty(ma.w3)
 if isempty(ma.w4)
  gradg=zeros(length(ma.w13c))
 #elseif !isempty(ma.w13c) #certaines variables sont fixés
 elseif !isempty(ma.w4) #certaines variables sont fixés
  tmp=zeros(ma.nb_comp)
  #tmp[ma.w3]=Relaxation.dpsi(xf[ma.w3+ma.n+ma.nb_comp],ma.r,ma.s,ma.t).*gradf[ma.w3+ma.n]
  tmp[ma.w4]=Relaxation.dpsi(xf[ma.w4+ma.n],ma.r,ma.s,ma.t).*gradf[ma.w4+ma.nb_comp+ma.n]
  gradg=redd(ma,tmp,ma.w13c)
 else #ma.w13c est vide
  gradg=Float64[]
 end

 gradh=Array{Float64,1}
 #if isempty(ma.w2) && isempty(ma.w4)
 if isempty(ma.w3)
  gradh=zeros(length(ma.w24c))
 #elseif !isempty(ma.w24c)
 elseif !isempty(ma.w3)
  tmp=zeros(ma.nb_comp)
  #tmp[ma.w4]=Relaxation.dpsi(xf[ma.w4+ma.n],ma.r,ma.s,ma.t).*gradf[ma.w4+ma.nb_comp+ma.n]
  tmp[ma.w3]=Relaxation.dpsi(xf[ma.w3+ma.n+ma.nb_comp],ma.r,ma.s,ma.t).*gradf[ma.w3+ma.n]
  gradh=redd(ma,tmp,ma.w24c)
 else
  gradh=Float64[]
 end

 return vcat(gradf[1:ma.n],gradf[ma.w13c+ma.n]+gradg,gradf[ma.w24c+ma.nb_comp+ma.n]+gradh)
end

"""
Evalue la matrice hessienne de la fonction objectif d'un MPCC actif
x est le vecteur réduit
"""
function hess(ma::MPCC_actif,x::Vector)

 #on calcul xf le vecteur complet
 xf=evalx(ma,x)

 #construction de la hessienne de taille (n+2nb_comp)^2

 #H=NLPModels.hess(ma.nlp,xf) #renvoi la triangulaire inférieure tril(H,-1)'
 #H=H+tril(H,-1)'

 H=ma.nlp.H(xf)
 H=H+tril(H,-1)'

 if ma.nb_comp>0
  return hess(ma,x,H)
 else
  return H
 end
end

"""
A partir de la hessienne complète (ou une approximation) : calcul la hessienne dans le sous-espace actif
"""
function hess(ma::MPCC_actif,x::Vector,H::Array{Float64,2})

 #on calcul xf le vecteur complet
 xf=evalx(ma,x)
 nred=length(x)
 nnb=ma.n+ma.nb_comp;nnbt=ma.n+2*ma.nb_comp;
 #construction du vecteur gradient de taille n+2nb_comp
 gradf=NLPModels.grad(ma.nlp,xf)

#la hessienne des variables du sous-espace (nredxnred)
 Hred=vcat(hcat(H[1:ma.n,1:ma.n],H[1:ma.n,ma.n+ma.w13c],H[1:ma.n,nnb+ma.w24c]),
             hcat(H[ma.n+ma.w13c,1:ma.n],H[ma.n+ma.w13c,ma.n+ma.w13c],H[ma.n+ma.w13c,nnb+ma.w24c]),
             hcat(H[nnb+ma.w24c,1:ma.n],H[nnb+ma.w24c,ma.n+ma.w13c],H[nnb+ma.w24c,ma.n+ma.nb_comp+ma.w24c]))

 if isempty(ma.w4)
  hessg=sparse(zeros(length(ma.w13c),nred))
 else
  hessg=sparse(zeros(length(ma.w13c),nred))

  w4r=zeros(Int64,length(ma.w4))
  for i=1:length(ma.w4)
   w4r[i]=findfirst(x->x==ma.w4[i],ma.w13c)
  end
  hessg=diagm(Relaxation.ddpsi(xf[ma.w4+ma.n],ma.r,ma.s,ma.t).*gradf[ma.w4+ma.nb_comp+ma.n])
  hessg+=diagm(Relaxation.dpsi(xf[ma.w4+ma.n],ma.r,ma.s,ma.t))*H[ma.w4+ma.nb_comp+ma.n,ma.w4+ma.nb_comp+ma.n]

  Hred[ma.n+w4r,ma.n+w4r]+=hessg
 end

 if isempty(ma.w3)
  hessh=sparse(zeros(length(ma.w24c),nred))
 else
  hessh=sparse(zeros(length(ma.w24c),nred))

  w3r=zeros(Int64,length(ma.w3))
  for i=1:length(ma.w3)
   w3r[i]=findfirst(x->x==ma.w3[i],ma.w24c)
  end
  hessh=diagm(Relaxation.ddpsi(xf[ma.w3+ma.nb_comp+ma.n],ma.r,ma.s,ma.t).*gradf[ma.w3+ma.n])
  hessh+=diagm(Relaxation.dpsi(xf[ma.w3+ma.nb_comp+ma.n],ma.r,ma.s,ma.t))*H[ma.w3+ma.n,ma.w3+ma.n]

  Hred[ma.n+length(ma.w13c)+w3r,ma.n+length(ma.w13c)+w3r]+=hessh
 end

 return Hred
end

"""
LSQComputationMultiplier(ma::MPCC_actif,x::Vector) :
ma MPCC_Actif
xj in n+2nb_comp
gradpen in n+2nb_comp

calcul la valeur des multiplicateurs de Lagrange pour la contrainte de complémentarité en utilisant moindre carré
"""
function LSQComputationMultiplier(ma::MPCC_actif,gradpen::Vector,xj::Vector)

 dg=Relaxation.dphi(xj[ma.n+1:ma.n+ma.nb_comp],xj[ma.n+ma.nb_comp+1:ma.n+2*ma.nb_comp],ma.r,ma.s,ma.t)
 gx=dg[1:ma.nb_comp];gy=dg[ma.nb_comp+1:2*ma.nb_comp]

 #matrices des contraintes actives : (lg,lh,lphi)'*A=b
 nw1=length(ma.w1)
 nw2=length(ma.w2)
 nwcomp=length(ma.wcomp)
 Dlg=-diagm(ones(nw1))
 Dlh=-diagm(ones(nw2))
 Dlphig=diagm(collect(gx)[ma.wcomp])
 Dlphih=diagm(collect(gy)[ma.wcomp])
 A=[hcat(Dlg,zeros(nw1,nw2+2*nwcomp));hcat(zeros(nw2,nw1),Dlh,zeros(nw2,2*nwcomp));hcat(zeros(nwcomp,nw1+nw2),Dlphig,Dlphih)] # nw1+nw2+nwcomp x nw1+nw2+2nwcomp
 #second membre
 b=-[gradpen[ma.w1+ma.n];gradpen[ma.wcomp+ma.n];gradpen[ma.w2+ma.n+ma.nb_comp];gradpen[ma.wcomp+ma.n+ma.nb_comp]] #nw1+nw2+2nwcomp
 #on calcule la solution par pseudo-inverse :
 l=pinv(A')*b

 lk=zeros(3*ma.nb_comp)
 lk[ma.w1]=l[1:nw1]
 lk[ma.nb_comp+ma.w2]=l[nw1+1:nw1+nw2]
 lk[2*ma.nb_comp+ma.wcomp]=l[nw1+nw2+1:nw1+nw2+nwcomp]

 return lk[1:ma.nb_comp],lk[ma.nb_comp+1:2*ma.nb_comp],lk[2*ma.nb_comp+1:3*ma.nb_comp]
end

"""
Définit la règle de relaxation de l'ensemble des contraintes
RelaxationRule(ma::ActifMPCCmod.MPCC_actif,xj::Vector,lg,lh,lphi,wmax)
xj          : de taille n+2nb_comp
(lg,lh,phi) : valeur multiplicateurs
wmax        : ensemble des contraintes qui viennent d'être ajouté

output : MPCC_actif (avec les ensembles de contraintes actives mis à jour
"""
function RelaxationRule(ma::ActifMPCCmod.MPCC_actif,xj::Vector,lg::Vector,lh::Vector,lphi::Vector,wmax::Array{Bool,2})

  copy_wmax=copy(wmax)

  # Relaxation de l'ensemble d'activation : désactive toutes les contraintes négatives
  ma.w[find(x -> x<0,[lg;lphi;lh;lphi])]=zeros(Bool,length(find(x -> x<0,[lg;lphi;lh;lphi])))
  # Règle d'anti-cyclage : on enlève pas une contrainte qui vient d'être ajouté.
  ma.w[find(x->x==1.0,copy_wmax)]=ones(Bool,length(find(x->x==1.0,copy_wmax)))

 return ActifMPCCmod.updatew(ma)
end

"""
Calcul le pas maximum que l'on peut prendre dans une direction d (par rapport à la contrainte de complémentarité relaxé)
d : direction réduit
xj : itéré réduit

output :
alpha : le pas maximum
w_save : l'ensemble des contraintes qui vont devenir actives si choisit alphamax
"""
function PasMaxComp(ma::MPCC_actif,x::Vector,d::Vector)

 #initialisation
 alpha=Inf #pas maximum que l'on peut prendre
 w_save=copy(ma.w) #double tableau des indices avec les contraintes activent en x+alpha*d

 #les indices où la première composante est libre
 for i in ma.w13c
  wr13=findfirst(x->x==i,ma.w13c) #l'indice relatif dans les variables libre
  iw13c=ma.n+wr13

  bloque=(i in ma.w2) && (i in ma.w4)
  if !(i in ma.w24c) && !bloque && d[iw13c]<0
   #on prend le plus petit entre x+alpha*dx>=-r et s+tTheta(x+alpha*dx-s)>=-r
   alpha11=(ma.nlp.meta.lvar[iw13c]-x[iw13c])/d[iw13c]
   alpha12=(Relaxation.invpsi(ma.nlp.meta.lvar[iw13c],ma.r,ma.s,ma.t)-x[iw13c])/d[iw13c]

   alphag=AlphaChoix(alpha,alpha11,alpha12)
   if alphag<=alpha
    alpha=alphag
    w_save=copy(ma.w)
   end

   #update of the active set
   if alpha11==alpha
    w_save[i,1]=true
   end
   if alpha12==alpha
    w_save[i+ma.nb_comp,2]=true
    w_save[i,2]=true
   end

  elseif bloque
   alpha=0.0
  end
 end

 #c'est un copie coller d'au dessus => exporter dans une fonction
 #les indices où la deuxième composante est libre
 for i in ma.w24c
  wr24=findfirst(x->x==i,ma.w24c) #l'indice relatif dans les variables libre
  iw24c=ma.n+length(ma.w13c)+wr24

  bloque=(i in ma.w1) && (i in ma.w3)
  if !(i in ma.w13c) && !bloque && d[iw24c]<0
   #on prend le plus petit entre y+alpha*dy>=-r et s+tTheta(y+alpha*dy-s)>=-r
   alpha21=(ma.nlp.meta.lvar[iw24c]-x[iw24c])/d[iw24c]
   alpha22=(Relaxation.invpsi(ma.nlp.meta.lvar[iw24c],ma.r,ma.s,ma.t)-x[iw24c])/d[iw24c]

   alphah=AlphaChoix(alpha,alpha21,alpha22)
   if alphah<=alpha
    alpha=alphah
    w_save=copy(ma.w)
   end

   #on met à jour les contraintes
   if alpha21==alpha
    w_save[i,2]=true
   end
   if alpha22==alpha
    w_save[i+ma.nb_comp,1]=true
    w_save[i,1]=true
   end
  elseif bloque
   alpha=0.0
  end
 end

 #enfin les indices où les deux sont libres
 for i in ma.wc
  #yG-psi(yH)=0 ou yH-psi(yG)=0
  wr1=findfirst(x->x==i,ma.w13c) #l'indice relatif dans les variables libre
  wr2=findfirst(x->x==i,ma.w24c) #l'indice relatif dans les variables libre
  iwr1=wr1+ma.n;iwr2=wr2+length(ma.w13c)+ma.n;

  #alphac=Relaxation.AlphaThetaMax(x[i+ma.n],d[i+ma.n],x[i+length(ma.w13c)+ma.n],d[i+length(ma.w13c)+ma.n],ma.r,ma.s,ma.t)
  alphac=Relaxation.AlphaThetaMax(x[iwr1],d[iwr1],x[iwr2],d[iwr2],ma.r,ma.s,ma.t)
  #yG-tb=0
  #alphac11=d[i+ma.n]<0 ? (ma.nlp.meta.lvar[ma.n+i]-x[i+ma.n])/d[i+ma.n] : Inf
  alphac11=d[iwr1]<0 ? (ma.nlp.meta.lvar[iwr1]-x[iwr1])/d[iwr1] : Inf
  #yH-tb=0
  #alphac21=d[i+length(ma.w13c)+ma.n]<0 ? (ma.nlp.meta.lvar[ma.n+length(ma.w13c)+i]-x[i+length(ma.w13c)+ma.n])/d[i+length(ma.w13c)+ma.n] : Inf
  alphac21=d[iwr2]<0 ? (ma.nlp.meta.lvar[iwr2]-x[iwr2])/d[iwr2] : Inf

  alphagh=AlphaChoix(alpha,alphac[1],alphac[2],alphac11,alphac21)

  if alphagh<=alpha
   alpha=alphagh
   w_save=copy(ma.w)
  end

   if alphac[1]==alpha
    w_save[i+ma.nb_comp,1]=true
   end
   if alphac[2]==alpha
    w_save[i+ma.nb_comp,2]=true
   end
   if alphac11==alpha
    w_save[i,1]=true
   end
   if alphac21==alpha
    w_save[i,2]=true
   end

 end #fin boucle for ma.wc

 return alpha,w_save,Array(w_save .& .!ma.w)
end

"""
Met à jour alpha si :
1) il est plus petit
2) il est non-nul
"""
function AlphaChoix(alpha::Float64,alpha1::Float64,alpha2::Float64)
 return AlphaChoix(alpha,alpha1,alpha2,0.0,0.0)
end

function AlphaChoix(alpha::Float64,alpha1::Float64,alpha2::Float64,alpha3::Float64,alpha4::Float64)
 prec=eps(Float64)
 a=alpha1,alpha2,alpha3,alpha4
 a=a[find(x->x>=prec,collect(a))]
 if isempty(a)
  a=max(alpha1,alpha2,alpha3,alpha4)
 end

 return min(minimum(a),alpha)
end

"""
Calcul le pas maximum que l'on peut prendre dans une direction d par rapport aux contraintes de bornes sur x.
d : direction réduit
xj : itéré réduit

output :
alpha : le pas maximum
w_save : l'ensemble des contraintes qui vont devenir actives si choisit alphamax
"""
function PasMaxBound(ma::MPCC_actif,x::Vector,d::Vector)
 if max(l-x)>0 || max(x-u)>0
  println("Error PasMaxBound: infeasible x")
 end
 #ma.nlp.meta.lvar[1:ma.n]
 #ma.nlp.meta.uvar[1:ma.n]
 return #TO DO : nécessite d'ajouter les x dans l'ensemble des contraintes actives
end

"""
Calcul le pas maximum que l'on peut prendre dans une direction d
d : direction réduit
xj : itéré réduit

output :
alpha : le pas maximum
w_save : l'ensemble des contraintes qui vont devenir actives si choisit alphamax
"""
function PasMax(ma::MPCC_actif,x::Vector,d::Vector)

 if ma.nb_comp>0
  #on récupère les infos sur la contrainte de complémentarité
  alpha,w_save,w_new=PasMaxComp(ma,x,d)
 else
  alpha=Inf;w_save=zeros(Bool,0,0);w_new=ma.w;
 end
 #alpha,w_save,w_new=PasMaxBound(ma,x,d)

 if alpha<0.0
  println("PasMax error: pas maximum négatif.")
  return
 end

 return alpha,w_save,w_new
end

#end of module
end


text:
"""
Runs the forward model for a given theta:

hs1d, taus1d, ivs1d, hs2d, taus2d, ivs2d, taus1d_l, vol_ratio, gb, bdot1d, dhdt1d = run_theta(theta, gl, gb, pm, pn, pp);
"""
function run_theta(theta, gl, gb, pm, pn, pp)
    pmm = MPara(pm, todict(theta))

    # this was missing:
    dhdt1d = [dhdtfn(b, pmm.dhdt_cutoff, pmm.dhdt_base) for b in gb.bands]
    bdot1d = map_onto_bands(gb, gl.bdot, mean, FILL)

    @time hs1d, taus1d, ivs1d, hs2d, taus2d, ivs2d, taus1d_l, vol_ratio, gb, bdot1d, dhdt1d =
        fwdm(gb, gl, pp, pmm, pn, bdot1d, dhdt1d)
    return hs1d, taus1d, ivs1d, hs2d, taus2d, ivs2d, taus1d_l, vol_ratio, gb, bdot1d, dhdt1d
end

plan = """

- MPara gives the parameters needed.
- gl+pm give an "instance" of a glacier to run the model over.
  - mostly pm would have information on how to change some of the 1D input variables
  - this would need some checking whether gb needs updating

"""

"""
Returns a dict of SHAs and branches.  Errors if one of the repos is dirty,
unless error_on_dirty=false.
"""
function get_git_shas(;error_on_dirty=true)
    LG = Base.LibGit2
    out = Dict{String,Any}()

    for rep in ["KissMCMC", "VAWTools", "BITEModel"]
        gr = LG.GitRepo(Pkg.dir(rep))
        error_on_dirty && LG.isdirty(gr) && error("KissMCMC is dirty.  Either commit or force with error_on_dirty=false")
        out[rep] = [string(LG.head_oid(gr)), string(LG.headname(gr))]
    end
    return out
end

# # patch JLD: Not needed for JLD2
# JLD.readas(::Val{:FunctionSerializer}) = (x...) -> error("Storing functions with JLD is not supported")
# @suppress JLD.writeas(fn::Function) = Val{:FunctionSerializer}()

"""
Saves an MCMC model run.
"""
function savemcmc(thetas, blobs, gb::Bands,
                  pp::Phys, pm::MPara, pn::Num, pmcmc, #::MCMCNum,
                  theta0, run::Symbol, sigmas::Tuple;
                  flname="", force_overwrite=false, error_on_dirty=true)
    gl = gb.gl
    sh, siv1d, siv2d = sigmas
    sh = round(Int,sh)
    siv2d = round(Int,siv2d)
    if flname==""
        region = getregion(gl)
        lab = getname(gl)
        flname = "output/$region-$lab-$run-$sh-$siv2d.jld"
    end
    if ~force_overwrite && isfile(flname)
        error("File $flname exists!  Aborting.  Use 'force_overwrite=true'.")
    end
    shas = get_git_shas(error_on_dirty=error_on_dirty)

    file = JLD.jldopen(flname, "w")
    try
        file["thetas"] = thetas
        file["blobs"] = blobs
        file["gb"] = gb
        file["pp"] = pp
        file["pm"] = pm
        file["pn"] = pn
        file["pmcmc"] = pmcmc
        file["theta0"] = theta0
        file["run"] = run
        file["sigmas"] = sigmas
        file["git_shas"] = shas
        file["timestamp"] = Dates.now()
    finally
        close(file)
    end
    return flname
end

"""
Load results.

better use @load

gl, gb, pp, pm, pn, res, git_shas, timestamp = loadmcmc(flname)

thetas_e, accept_ratio_e, blobs_e = squash_chains(res... , drop_low_accept_ratio=true, blob_reduce! = blob_reduce!)
mc_vol, mc_vol_above, mc_hs1d, mc_ivs1d, mh2d, sh2d, miv2d, siv2d, pms = unpack_blob(blobs_e, size(gl.dem.v))
"""
function loadmcmc(flname)
    @unpack thetas, blobs, gb, pp, pm, pn, pmcmc,
    theta0, run, sigmas, git_shas, timestamp = JLD.load(flname)
    return thetas, blobs, gb, pp, pm, pn, pmcmc,
    theta0, run, sigmas, git_shas, timestamp
end

ext: jl


text:
function get_coefficients(obj::ShortPeriodTerms, arg0::AbsoluteDate, arg1::Set)
    return jcall(obj, "getCoefficients", Map, (AbsoluteDate, Set), arg0, arg1)
end

function get_coefficients_key_prefix(obj::ShortPeriodTerms)
    return jcall(obj, "getCoefficientsKeyPrefix", JString, ())
end

function value(obj::ShortPeriodTerms, arg0::Orbit)
    return jcall(obj, "value", Vector{jdouble}, (Orbit,), arg0)
end


text:
module LiveServer

# from stdlib
using Sockets, Pkg
# the only dependency (see the patch in http_patch.jl)
using HTTP

export serve, servedocs

#
# Environment variables
#

"""Script to insert on a page for live reload, see `client.html`."""
const BROWSER_RELOAD_SCRIPT = read(joinpath(@__DIR__, "client.html"), String)

"""Whether to display messages while serving or not, see [`verbose`](@ref)."""
const VERBOSE = Ref{Bool}(false)

"""The folder to watch, either the current one or a specified one."""
const CONTENT_DIR = Ref{String}("")

"""List of files being tracked with WebSocket connections."""
const WS_VIEWERS = Dict{String,Vector{HTTP.WebSockets.WebSocket}}()

"""Keep track of whether an interruption happened while processing a websocket."""
const WS_INTERRUPT = Base.Ref{Bool}(false)

# to fix lots of false error messages from HTTP
# https://github.com/JuliaWeb/HTTP.jl/pull/546
# we do HTTP.Stream{HTTP.Messages.Request,S} instead of just HTTP.Stream to prevent the Julia warning about incremental compilation
# This hack was kindly suggested by Fons van der Plas, the author of Pluto.jl see
# https://github.com/fonsp/Pluto.jl/commit/34d41e63138ee6dad178cd9916d4721441eaf710
function HTTP.closebody(http::HTTP.Stream{HTTP.Messages.Request,S}) where S <: IO
    http.writechunked || return
    http.writechunked = false
    try; write(http.stream, "0\r\n\r\n"); catch; end
    return
end

#
# Functions
#

include("mimetypes.jl")
include("file_watching.jl")
include("server.jl")

include("utils.jl")

end # module


text:
"""
    gaussian_stats([::Type{T},] num_features::Int[, weight])

This function creates an OnlineStats.KahanVariance object for tracking the mean and variance for a vector.
Any OnlineStats.Weight can be used. The default is OnlineStats.EqualWeight and OnlineStats.ExponentialWeight
if an integer or float is given as the weight.

See also: [`extrema_stats`](@ref), [`LinearNormalization`](@ref)

# Examples
```jldoctest
julia> stats = gaussian_stats(Float32, 2, 1e-4)
Group
├─ KahanVariance: n=0 | value=1.0
└─ KahanVariance: n=0 | value=1.0

julia> fit!(stats, [1.0, 2.0])
Group
├─ KahanVariance: n=1 | value=1.0
└─ KahanVariance: n=1 | value=1.0

```
"""
function gaussian_stats end

function gaussian_stats(::Type{T}, num_features::Int, weight) where {T}
    return num_features*KahanVariance(T=T, weight=weight)
end

function gaussian_stats(::Type{T}, num_features::Int) where {T}
    return num_features*KahanVariance(T)
end

function gaussian_stats(::Type{T}, num_features::Int, weight::TW) where {T, TW<:Real}
    return gaussian_stats(T, num_features, ExponentialWeight(weight))
end

function gaussian_stats(num_features::Int, weight)
    return gaussian_stats(Float64, num_features, weight)
end

function gaussian_stats(num_features::Int)
    return gaussian_stats(Float64, num_features)
end

"""
    extrema_stats([::Type{T},] num_features::Int)

This function creates an OnlineStats.KahanVariance object for tracking the mean and variance for a vector.
Any OnlineStats.Weight can be used. The default is OnlineStats.EqualWeight and OnlineStats.ExponentialWeight
if an integer or float is given as the weight.

See also: [`extrema_stats`](@ref), [`LinearNormalization`](@ref)

# Examples
```jldoctest
julia> stats = gaussian_stats(Float32, 2, 1e-4)
Group
├─ KahanVariance: n=0 | value=1.0
└─ KahanVariance: n=0 | value=1.0

julia> fit!(stats, [1.0, 2.0])
Group
├─ KahanVariance: n=1 | value=1.0
└─ KahanVariance: n=1 | value=1.0

```
"""
function extrema_stats(::Type{T}, num_features::Int) where {T}
    return num_features*Extrema(T)
end

function extrema_stats(num_feature::Int)
    return num_feature*Extrema(Float64)
end


text:
mutable struct childSupportIncome:>
    nonfixedIncome
end

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
using MLBase
# from Taddy (2015, distrom):
# CV selection rules: both CV1se, which chooses the largest λt with mean
# OOS deviance no more than one standard error away from minimum, and CVmin,
# which chooses λt at lowest mean OOS deviance.

function CVmin(oosdevs)
    cvmeans = vec(mean(oosdevs,dims=2))
    segCVmin = argmin(cvmeans)
end

function CV1se(oosdevs)
    nλ,nfolds = size(oosdevs)
    cvmeans = vec(mean(oosdevs,dims=2))
    (mincvmean,segCVmin) = findmin(cvmeans)
    mincvstds = std(view(oosdevs, segCVmin, :)) / sqrt(nfolds-1)
    mincvmean_plus_mincvstds = mincvmean + mincvstds
    for s=1:nλ
        cv1se = mincvmean_plus_mincvstds - cvmeans[s]
        if cv1se >= 0
            return s
        end
    end
    error("should have found the cv1se by now")
end

# convenience function to use the same data as in original path
function cross_validate_path(path::RegularizationPath;     # fitted path
                        gen=Kfold(length(path.m.rr.y),10), # folds generator (see MLBase)
                        select=:CVmin,                     # :CVmin or :CV1se
                        fitargs...)
    m = path.m
    y = m.rr.y
    offset = m.rr.offset
    Xstandardized = m.pp.X
    cross_validate_path(path,Xstandardized,y;gen=gen,select=select,offset=offset,
        standardize=false,fitargs...)
end

pathtype(::LassoPath) = LassoPath
pathtype(::GammaLassoPath) = GammaLassoPath

function cross_validate_path(path::RegularizationPath,    # fitted path
                       X::AbstractMatrix{T}, y::V;        # potentially new data
                       gen=Kfold(length(y),10),           # folds generator (see MLBase)
                       select=:CVmin,                     # :CVmin or :CV1se
                       offset::FPVector=T[],
                       fitargs...) where {T<:AbstractFloat,V<:FPVector}
    @extractfields path m λ
    n,p = size(X)
    @assert n == length(y) "size(X,1) != length(y)"

    nfolds = length(gen)
    nλ = length(λ)
    d = distfun(path)
    l = linkfun(path)

    # valid offset given?
    if length(m.rr.offset) > 0
        length(offset) == n ||
            throw(ArgumentError("fit with offset, so `offset` kw arg must be an offset of length $n"))
    else
        length(offset) > 0 && throw(ArgumentError("fit without offset, so value of `offset` kw arg does not make sense"))
    end

    # EQUAL WEIGHTS ONLY!
    wts = ones(T, n)

    # results array
    oosdevs = zeros(T,nλ,nfolds)

    for (f, train_inds) in enumerate(gen)
        test_inds = setdiff(1:n, train_inds)
        nis = length(test_inds)

        if length(offset) > 0
            foldoffset = offset[train_inds]
        else
            foldoffset = offset
        end

        # fit model to train_inds
        foldpath = fit(pathtype(path),X[train_inds,:],y[train_inds],d,l;λ=λ,wts=wts[train_inds],offset=foldoffset,fitargs...)

        if length(offset) > 0
            foldoffset = offset[test_inds]
        else
            foldoffset = offset
        end

        # calculate etas for each obs x segment
        μ = predict(foldpath, X[test_inds,:]; offset=foldoffset, select=:all)

        # calculate deviations on test sets efficiently (not much mem)
        for s=1:nλ
            # deviance of segment s (cummulator for sum of obs deviances)
            devs = zero(T)

            for ip=1:nis
                # get test obs
                yi = y[test_inds[ip]]

                # deviance of a single observation i in segment s
                devs += devresid(d, yi, μ[ip,s])
            end

            # store result
            oosdevs[s,f] = devs/nis
        end
    end

    CVfun = eval(select)
    CVfun(oosdevs)
end


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module GeneralizedBinomial
using Nemo, PrimesIterator, NumberTheory, Products, Triangles

export ModuleGeneralizedBinomial
export Binomial, Pascal, T007318

"""
 P. Luschny, [Generalized Binomial](http://oeis.org/wiki/User:Peter_Luschny/ExtensionsOfTheBinomial), OEIS Wiki.

* Binomial, Pascal, T007318
"""
const ModuleGeneralizedBinomial = ""

"""
The classical binomial coefficients defined for ``n≥0`` and ``0≤k≤n`` (a.k.a. Pascal's triangle).
"""
function Pascal(n::Int, k::Int)

    (k == 0 || k == n) && return 1
    if k > div(n, 2) k = n - k end

    nk = n - k
    factors = fmpz[]
    rootN = isqrt(n)

    # Make use of Kummer's theorem.
    for prime in Primes(2, n)

        if prime > nk
            push!(factors, prime)
            continue
        end

        prime > div(n, 2) && continue

        if prime > rootN
            (n % prime < k % prime) && push!(factors, prime)
            continue
        end

        r, N, K, p = 0, n, k, 1

        while N > 0
            r = N % prime < (K % prime + r) ? 1 : 0
            if r == 1 p *= prime end
            N = div(N, prime)
            K = div(K, prime)
        end

        p > 1 && push!(factors, p)
    end

    ∏(factors) end

"""
Pascal's triangle.
"""
function T007318(n::Int)
    T = zeros(QQ, div(n * (n + 1), 2))
    j = 1
    for m in 0:n - 1, k in 0:m
        T[j] = binom(m, k)
        j += 1
    end
    T
end

# See the discussion on
# [Extensions of the Binomial](http://oeis.org/wiki/User:Peter_Luschny/ExtensionsOfTheBinomial).
"""
Return the extended binomial coefficients defined for all ``n ∈ Z`` and ``k ∈ Z``. Behaves like the binomial function in Maple and Mathematica.

``\\binom{n}{k} = \\lim \\limits_{x \\rightarrow 1}(Γ(n + x) / (Γ(k + x) Γ(n - k + x))).``

"""
function Binomial(n::Int, k::Int)
    0 ≤ k ≤ n  && return binom(n, k)
    k ≤ n <  0 && return binom(-k - 1, n - k) * (-1)^(n - k)
    n <  0 ≤ k && return binom(-n + k - 1, k) * (-1)^k
    ZZ(0)
end

#START-TEST-########################################################

using Test
function test()

    @testset "Binomial" begin
        for n in 0:10, k in 0:n
            @test Binomial(n, k) == div(fac(n), (fac(n - k) * fac(k)))
            @test Binomial(n, k) == Pascal(n, k)
        end
    end
end

function demo()
    ShowAsΔ(T007318(8))
    println()

    for n in -10:10
        println([Binomial(n, k) for k in -10:10])
    end
end

"""
for n in 0:10000 Binomial(2*n,n) end
    0.504729 seconds (10.00 k allocations: 156.266 KiB)
for n in -100:100, k in -100:100 Binomial(n,k) end
    0.008669 seconds (55.55 k allocations: 867.984 KiB)
for k in -10000:10000 Binomial(-5,k) end
    0.005378 seconds (40.00 k allocations: 624.969 KiB)
"""
function perf()
    @time (for n in 0:10000 Binomial(2n, n) end)
    @time (for n in -100:100, k in -100:100 Binomial(n, k) end)
    @time (for k in -10000:10000 Binomial(-5, k) end)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module


text:
using Compat.readline

# find a good port
conn = listenany(2228)
close(conn[2])
port = conn[1]

server = @async lintserver(port)
sleep(1) #let server start

function lintbyserver(socket,str)
    println(socket, "none")
    println(socket, sizeof(str)) # bytes of code
    println(socket, str) # code
end

function readfromserver_old(socket)
    response = ""
    line = ""
    while line != "\n"
        response *= line
        line = readline(socket)
    end
    return response
end

function readfromserver_new(socket)
    response = ""
    line = ""
    while isopen(socket)
        response *= line
        line = readline(socket)
    end
    return response
end

@testset "lintserver() tests" begin
    conn = connect(port)
    write(conn, "empty\n")
    write(conn, "1\n")
    write(conn, "\n")

    @test readline(conn) == ""

    conn = connect(port)
    write(conn, "undeclared_symbol\n")
    write(conn, "4\n")
    write(conn, "bad\n")

    @test occursin("use of undeclared symbol", readline(conn))
    @test readline(conn) == ""
end

@testset "Testing the lintserver addition" begin
    str = """
          test = "Hello" + "World"
          """
    socket = connect(port)
    lintbyserver(socket, str)
    response = String[]
    line = "."
    while !isempty(line)
        line = readline(socket)
        push!(response, line)
    end

    @test "none:1 E422 : string uses * to concatenate" in response

    socket = connect(port)
    lintbyserver(socket, str)
    response = String[]
    while isopen(socket)
        line = readline(socket)
        push!(response, line)
    end

    @test "none:1 E422 : string uses * to concatenate" in response
end

@testset "Testing lintserver() with named pipe and JSON format" begin
    if is_windows()
        pipe_lm = "\\\\.\\pipe\\testsocket"
        pipe_slv1 = "\\\\.\\pipe\\testsocket2"
        pipe_vscode = "\\\\.\\pipe\\testsocket3"
        pipe_slv2 = "\\\\.\\pipe\\testsocket4"
    else
        pipe_lm = tempname()
        pipe_slv1 = tempname()
        pipe_vscode = tempname()
        pipe_slv2 = tempname()
    end

    function writeandreadserver(pipe,json_input)
        conn = connect(pipe)
        JSON.print(conn, json_input)
        JSON.parse(conn)
    end

    server_lm = @async lintserver(pipe_lm,"lint-message")
    server_slv1 = @async lintserver(pipe_slv1,"standard-linter-v1")
    server_vscode = @async lintserver(pipe_vscode,"vscode")
    server_slv2 = @async lintserver(pipe_slv2,"standard-linter-v2")
    sleep(1)

    json_input1 = Dict("file" => "none", "code_str" => "something")
    json_input2 = Dict("file" => "none", "code_str" => "pi=3")
    json_input3 = Dict("file" => "none", "code_str" => "function a(b)\nend")

    results_array = writeandreadserver(pipe_lm, json_input1)
    @test_broken results_array[1]["line"] == 1
    @test results_array[1]["message"] == "use of undeclared symbol"
    @test_broken results_array[1]["file"] == "none"
    @test results_array[1]["code"] == "E321"

    results_array = writeandreadserver(pipe_slv1, json_input1)
    @test results_array[1]["text"] == "E321 something: use of undeclared symbol"
    @test results_array[1]["filePath"] == "none"
    @test results_array[1]["range"] == Array[[0, 0], [0, 80]]
    @test results_array[1]["type"] == "error"

    results_array = writeandreadserver(pipe_slv1, json_input2)
    @test startswith(results_array[1]["text"], "I343 pi: global variable")
    @test results_array[1]["filePath"] == "none"
    @test results_array[1]["range"] == Array[[0, 0], [0, 80]]
    @test results_array[1]["type"] == "info"

    results_array = writeandreadserver(pipe_slv1, json_input3)
    @test startswith(results_array[1]["text"], "I340 b: unused local variable")
    @test results_array[1]["filePath"] == "none"
    @test_broken results_array[1]["range"] == Array[[0, 0], [0, 80]]
    @test results_array[1]["type"] == "info"

    results_array = writeandreadserver(pipe_vscode, json_input1)
    @test results_array[1]["message"] == "something: use of undeclared symbol"
    @test results_array[1]["filePath"] == "none"
    @test results_array[1]["range"] == Array[[0, 0], [0, 80]]
    @test results_array[1]["code"] == "E321"
    @test results_array[1]["severity"] == 1
    @test results_array[1]["source"] == "Lint.jl"

    results_array = writeandreadserver(pipe_slv2, json_input1)
    @test results_array[1]["description"] == "E321"
    @test results_array[1]["location"]["file"] == "none"
    @test results_array[1]["location"]["position"] == Array[[0, 0], [0, 80]]
    @test results_array[1]["severity"] == "error"
    @test results_array[1]["excerpt"] == "something: use of undeclared symbol"

    json_input4 = Dict("file" => "none", "code_str" => "function a(b)\nend",
                       "ignore_info" => true)
    results_array = writeandreadserver(pipe_lm, json_input4)
    @test isempty(results_array)

    json_input5 = Dict("file" => "none", "code_str" => "pi = 1",
                       "ignore_warnings" => true)
    results_array = writeandreadserver(pipe_lm, json_input5)
    @test_broken isempty(results_array)

    json_input6 = Dict("file" => "none",
                       "code_str" => "pi = 1\nfunction a(b)\nend",
                       "ignore_codes" => ["I340","I343"])
    results_array = writeandreadserver(pipe_lm, json_input6)
    @test isempty(results_array)

    json_input7 = Dict("file" => "none",
                       "code_str" => "pi = 1\nfunction a(b)\nend",
                       "show_code" => false)
    results_array = writeandreadserver(pipe_slv1, json_input7)
    @test startswith(results_array[1]["text"], "pi: global variable")
    @test results_array[1]["filePath"] == "none"
    @test_broken results_array[1]["range"] == Array[[0, 0], [0, 80]]
    @test_broken results_array[1]["type"] == "warning"

    json_input8 = Dict("file" => "none",
                       "code_str" => "pi = 1\nfunction a(b)\nend",
                       "show_code" => true)
    results_array = writeandreadserver(pipe_slv1, json_input8)
    @test startswith(results_array[1]["text"], "I343 pi: global variable")
    @test results_array[1]["filePath"] == "none"
    @test_broken results_array[1]["range"] == Array[[0, 0], [0, 80]]
    @test results_array[1]["type"] == "info"
end


text:
module day03
using CircularArrays

include("inputs.jl")
ex = inputs.ex_d3()
input = inputs.input_d3()

# Part 1
function parsemap(input::String)
    map = split(input) .|> # array of strings "..##......."
    collect .|> # array of array of chars ['.', '.', '#', '#', ...]
    x->(x .=='#') |> # array of BitArray [0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0]
    x->hcat(x...)    # array of Bools 0  0  1  1  0  0  0  0  0  0  0
    map = vcat(map...) |>
    CircularArray
end

function count_trees(input, slope=(3,1))
    map = parsemap(input)
    count = 0
    x, y = 1, 1
    while y<=size(map)[1]
        count += map[y, x] ? 1 : 0
        x += slope[1]
        y += slope[2]
    end
    return count
end

@assert count_trees(ex) == 7
@info "Day 3, Part 1 answer: $(count_trees(input))"

# Part 2
function assess_slopes(input)
    slopes = [(1,1), (3,1), (5,1), (7,1), (1,2)]
    counts = [count_trees(input, slope) for slope in slopes]
    return reduce(*, counts)
end

@assert assess_slopes(ex) == 336
@info "Day 3, Part 2 answer: $(assess_slopes(input))"

end


text:
using StatsBase
using YaoArrayRegister
export UnitaryChannel

# TODO: replace this with uweight
# when StatsBase 0.33 is out
_uweights(n) = weights(ones(n))

# NOTE: this can have a better support in YaoIR directly with StatsBase
# as an compiled instruction, but store the weights is the best solution
# here

"""
    UnitaryChannel(operators[, weights])

Create a unitary channel, optionally weighted from an list of weights.
The unitary channel is defined as below in Kraus representation

```math
ϕ(ρ) = \\sum_i U_i ρ U_i^†
```

!!! note
    Unitary channel will only normalize the weights when calculating the matrix form,
    thus you should be careful when you need this condition for other purpose.

!!! note
    when applying a `UnitaryChannel` on the register, a unitary will be sampled
    uniformly or optionally from given weights, then this unitary will be applied
    to the register.

# Example

```jldoctest; setup=:(using YaoBlocks, YaoArrayRegister)
julia> UnitaryChannel([X, Y, Z])
nqudits: 1
unitary_channel
├─ [1.0] X
├─ [1.0] Y
└─ [1.0] Z
```

Or with weights

```jldoctest; setup=:(using YaoBlocks, YaoArrayRegister)
julia> UnitaryChannel([X, Y, Z], [0.1, 0.2, 0.7])
nqudits: 1
unitary_channel
├─ [0.1] X
├─ [0.2] Y
└─ [0.7] Z
```
"""
struct UnitaryChannel{W<:AbstractWeights} <: CompositeBlock{2}
    n::Int
    operators::Vector{AbstractBlock{2}}
    weights::W

    function UnitaryChannel(operators::Vector)
        w = _uweights(length(operators))
        n = _check_block_sizes(operators)
        new{typeof(w)}(n, operators, w)
    end

    function UnitaryChannel(operators::Vector, w::AbstractWeights)
        n = _check_block_sizes(operators)
        new{typeof(w)}(n, operators, w)
    end

    function UnitaryChannel(operators::Vector, w::AbstractVector)
        w = weights(w)
        n = _check_block_sizes(operators)
        new{typeof(w)}(n, operators, w)
    end
end

UnitaryChannel(it, weights) = UnitaryChannel(collect(it), weights)
UnitaryChannel(it) = UnitaryChannel(collect(it))
nqudits(uc::UnitaryChannel) = uc.n

function _apply!(r::AbstractRegister, x::UnitaryChannel)
    _apply!(r, sample(x.operators, x.weights))
end

function mat(::Type{T}, x::UnitaryChannel) where {T}
    U = sample(x.operators, x.weights)
    return mat(T, U)
end

subblocks(x::UnitaryChannel) = x.operators
chsubblocks(x::UnitaryChannel, it) = UnitaryChannel(collect(it), x.weights)
occupied_locs(x::UnitaryChannel) = union(occupied_locs.(x.operators)...)

function cache_key(x::UnitaryChannel)
    key = hash(x.weights)
    for each in x.operators
        key = hash(each, key)
    end
    return key
end

function Base.:(==)(lhs::UnitaryChannel, rhs::UnitaryChannel)
    return (lhs.n == rhs.n) && (lhs.weights == rhs.weights) && (lhs.operators == rhs.operators)
end

Base.adjoint(x::UnitaryChannel) = UnitaryChannel(adjoint.(x.operators), x.weights)


text:
function _stripcomments(txt::AbstractString, dlm::Tuple{String,String})
    "Strips first nest of block comments"

    dlml, dlmr = dlm
    indx = searchindex(txt, dlml)
    if indx > 0
        out = IOBuffer()
        write(out, txt[1:indx-1])
        txt = txt[indx+length(dlml):end]
        txt = _stripcomments(txt, dlm)
        indx = searchindex(txt, dlmr)
        @assert(indx > 0, "cannot find a closer delimiter \"$dlmr\" in $txt")
        write(out, txt[indx+length(dlmr):end])
    else
        out = txt
    end
    return String(out)
end

function stripcomments(txt::AbstractString, dlm::Tuple{String,String}=("/*", "*/"))
    "Strips nests of block comments"

    dlml, dlmr = dlm
    while contains(txt, dlml)
        txt = _stripcomments(txt, dlm)
    end

    return txt
end

function main()
    println("\nNON-NESTED BLOCK COMMENT EXAMPLE:")
    smpl = """
/**
* Some comments
* longer comments here that we can parse.
*
* Rahoo
*/
function subroutine() {
a = /* inline comment */ b + c ;
}
/*/ <-- tricky comments */

/**
* Another comment.
*/
function something() {
}
"""
    println(stripcomments(smpl))

    println("\nNESTED BLOCK COMMENT EXAMPLE:")
    smpl = """
/**
* Some comments
* longer comments here that we can parse.
*
* Rahoo
*//*
function subroutine() {
a = /* inline comment */ b + c ;
}
/*/ <-- tricky comments */
*/
/**
* Another comment.
*/
function something() {
}
"""
    println(stripcomments(smpl))
end

main()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# A version of evolution with selection for robustness that evolves population until all individuals are at optimal goal fitness,
#    then evolves with selection for mutational robustness for a fixed number of generations.
# A limitation is that fitness is not maintained during mutation for robustenss.
export run_robust_evolve, robust_evolve
using DataFrames
using Statistics
using CSV
using Distributed
using Main.CGP

# To run julia:  julia -i -p 4 -L CGP.jl

# On 8/19/20 commented out the "var" fields of the dataframe

function run_robust_evolve( nreps::Int64, p::Parameters, popsize::Int64, max_pop_gens::Int64, indiv_steps::Int64, robust_steps::Int64, fit_steps::Int64,
      goallist::GoalList, tourn_size::Int64, csvfile::String )   # tourn_size==0 means proportional selection
  funcs = default_funcs(p.numinputs)
  test_MyInt(p.numinteriors)
  df_results = DataFrame[]
  df = DataFrame()
  df.generation = collect(1:2*robust_steps)
  df.numinputs = fill(p.numinputs,2*robust_steps)
  df.numoutputs = fill(p.numoutputs,2*robust_steps)
  df.numints = fill(p.numinteriors,2*robust_steps)
  df.levsback = fill(p.numlevelsback,2*robust_steps)
  df.popsize = fill(popsize,2*robust_steps)
  df.robust_sel = fill(false,2*robust_steps)
  df.mean_robfit = zeros(Float64,2*robust_steps)
  df.max_robfit = zeros(Float64,2*robust_steps)
  df.mean_nactive = zeros(Float64,2*robust_steps)
  #df.var_nactive = zeros(Float64,2*robust_steps)
  df.mean_complexity = zeros(Float64,2*robust_steps)
  #df.var_complexity = zeros(Float64,2*robust_steps)
  df.mean_degeneracy = zeros(Float64,2*robust_steps)
  #df.var_degeneracy = zeros(Float64,2*robust_steps)
  df_init = [ deepcopy(df) for _ = 1:nreps]
  df_results = pmap( ddf->robust_evolve(ddf,p, popsize, max_pop_gens, indiv_steps, robust_steps, fit_steps, goallist, tourn_size ), df_init)
  #df_results = map( ddf->robust_evolve(ddf,p, popsize, max_pop_gens, indiv_steps, robust_steps, fit_steps, goallist, tourn_size ), df_init)
  #println("df_results.mean_complexity: ",[df_results[i].mean_complexity for i = 1:nreps])
  #println("df_results.mean_complexity: ",[d.mean_complexity for d in df_results])
  df.robust_sel = df_results[1].robust_sel
  # Note that df.mean_robfit is a vector over reps of vectors over generations
  df.mean_robfit = mean([d.mean_robfit for d in df_results])
  df.max_robfit = mean([d.max_robfit for d in df_results])
  # Note that df.mean_nactive is a vector over reps of vectors over generations
  # the mean() function computes the mean over reps for each generation
  # the var() function computes the sample variance over reps for each generation
  df.mean_nactive = mean([d.mean_nactive for d in df_results])
  @assert df.mean_nactive == sum([d.mean_nactive for d in df_results])/nreps
  #df.var_nactive = var([d.mean_nactive for d in df_results])
  df.mean_complexity = mean([d.mean_complexity for d in df_results])
  #df.var_complexity = var([d.mean_complexity for d in df_results])
  #println("df.var_complexity: ",df.var_complexity)
  df.mean_degeneracy = mean([d.mean_degeneracy for d in df_results])
  #df.var_degeneracy = var([d.mean_degeneracy for d in df_results])
  open( csvfile, "w" ) do f
    println(f,"MyInt: ",Main.CGP.MyInt)
    println(f,"funcs: ",funcs)
    println(f,"nreps: ",nreps)
    println(f,"maxints_for_degen: ",Main.CGP.maxints_for_degen)
    println(f,"indiv_steps: ",indiv_steps)
    println(f,"robust_steps: ",robust_steps)
    println(f,"fit_steps: ",fit_steps)
    println(f,"popsize: ",popsize)
    println(f,"max_pop_gens: ",max_pop_gens)
    println(f,"ngoals: ",length(goallist))
    CSV.write( f, df, append=true, writeheader=true )
  end
  df
end

function robust_evolve( df::DataFrame, p::Parameters, popsize::Int64, max_pop_gens::Int64, indiv_steps::Int64, robust_steps::Int64, fit_steps::Int64,
    goallist::GoalList, tourn_size::Int64 )   # tourn_size==0 means proportional selection
  funcs = default_funcs(p.numinputs)
  maxints_for_degen = Main.CGP.maxints_for_degen
  for robust_sel = false:true
    pop = [ random_chromosome(p,funcs) for _ = 1:popsize ]
    gg = 1
    while gg < max_pop_gens
      #println("gen: ",gg,"  fitnesses: ",[pop[i].fitness for i =1:popsize] )
      for i = 1:popsize
        (pop[i],step,worse,same,better,output,matched_goals,matched_goals_list) = mut_evolve( pop[i], goallist, funcs, indiv_steps )
        #println("(gg,i): ",(gg,i),"  pop[i].fitness: ",pop[i].fitness)
      end
      fitness_vector = [ pop[i].fitness for i = 1:popsize ]
      if minimum(fitness_vector) == Float64(p.numoutputs)
        break
      end
      gg += 1
    end
    if gg == max_pop_gens
      error("generation_limit_exceeded")
    else
      println("population optimized at genertion ",gg)
      #println("fitvec: ", [ pop[i].fitness for i = 1:popsize ])
    end
    gr_increment = robust_sel*robust_steps
    #println("gr_increment: ",gr_increment)
    for gr = (gr_increment + 1):(gr_increment + robust_steps)
      df.robust_sel[gr] = robust_sel
      for i = 1:popsize
        if robust_sel
          pop[i].robustness = mutational_robustness( pop[i], funcs )
          robfit_vector = [ pop[i].robustness for i = 1:popsize ]
        else
          (pop[i],step,worse,same,better,output,matched_goals,matched_goals_list) = mut_evolve( pop[i], goallist, funcs, fit_steps, hamming_sel=robust_sel)
          robfit_vector = [ pop[i].fitness for i = 1:popsize ]
        end
      end
      #println("robfit before selection: ", [ round2(v) for v in robfit_vector] )
      if tourn_size == 0
        propsel!( pop, robfit_vector, maxfit=findmax(robfit_vector)[1] )
      else
        tournsel!( pop, robfit_vector, tourn_size )
      end
      #println("df.numints: ",df.numints)
      robfit_vector = robust_sel ? [ pop[i].robustness for i = 1:popsize ] : [ pop[i].fitness for i = 1:popsize ]
      nactive_vector = [ number_active( pop[i] ) for i = 1:popsize ]
      complexity_vector = df.numints[1] <= maxints_for_degen ? [ complexity5( pop[i] ) for i = 1:popsize ] : zeros(MyInt,popsiz)
      degeneracy_vector = df.numints[1] <= maxints_for_degen ? [ degeneracy( pop[i] ) for i = 1:popsize ] : zeros(MyInt,popsize)
      df.mean_robfit[gr] = mean(robfit_vector)
      df.max_robfit[gr] = maximum(robfit_vector)
      df.mean_nactive[gr] = mean(nactive_vector)
      df.mean_complexity[gr] = mean(complexity_vector)
      df.mean_degeneracy[gr] = mean(degeneracy_vector)
      #println("robfit after selection: ", [ round2(v) for v in robfit_vector] )
      #=
      if gr==7 || gr==14
        println("gr: ",gr,"  complexity: ",complexity_vector,"  mean compleixty: ",mean(complexity_vector))
      end
      =#
    end
    #robfit_vector = [ pop[i].robfit for i = 1:popsize ]
  end
  #println("robust_evolve finished")
  df
end


text:
import CUDA.CUDNN: cudnnPoolingForwardAD
using CUDA.CUDNN: cudnnPoolingBackward, handle
using AutoGrad: AutoGrad, @primitive1, value

@primitive1((cudnnPoolingForwardAD(x; poolingDesc, alpha, beta, xDesc, yDesc, y),
             _dy,_y),
            ((x,y,dy,dx) = (value(x),value(_y),value(_dy),similar(x));
             cudnnPoolingBackward(handle(), poolingDesc, alpha, yDesc, y, yDesc, dy, xDesc, x, beta, xDesc, dx);
             dx))


text:
immutable NoncentralBeta <: ContinuousUnivariateDistribution
    alpha::Float64
    beta::Float64
    ncp::Float64
    function NoncentralBeta(a::Real, b::Real, nc::Real)
    	a > 0.0 && b > 0.0 && nc >= 0.0 ||
            error("alpha and beta must be > 0 and ncp >= 0")
    	new(float64(a), float64(b), float64(nc))
    end
end

@_jl_dist_3p NoncentralBeta nbeta

insupport(::NoncentralBeta, x::Real) = zero(x) < x < one(x)
insupport(::Type{NoncentralBeta}, x::Real) = zero(x) < x < one(x)

function rand(d::NoncentralBeta)
    a = rand(NoncentralChisq(2.0 * d.alpha, d.ncp))
    b = rand(Chisq(2.0 * d.beta))
    a / (a+b)
end


text:
@testset "Variable Tests" begin
  v = Variable("x0", :ContinuousScalar)
  # NOTE: A timestamp substitution is done in these JSON payloads.
  @test JSON.json(v) == "{\"label\":\"x0\",\"dataEntry\":\"{}\",\"nstime\":\"0\",\"variableType\":\"IncrementalInference.ContinuousScalar\",\"dataEntryType\":\"{}\",\"ppeDict\":\"{}\",\"solverDataDict\":\"{\\\"default\\\":{\\\"vecval\\\":[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],\\\"dimval\\\":1,\\\"vecbw\\\":[0.0],\\\"dimbw\\\":1,\\\"BayesNetOutVertIDs\\\":[],\\\"dimIDs\\\":[0],\\\"dims\\\":1,\\\"eliminated\\\":false,\\\"BayesNetVertID\\\":\\\"_null\\\",\\\"separator\\\":[],\\\"variableType\\\":\\\"IncrementalInference.ContinuousScalar\\\",\\\"initialized\\\":false,\\\"infoPerCoord\\\":[0.0],\\\"ismargin\\\":false,\\\"dontmargin\\\":false,\\\"solveInProgress\\\":0,\\\"solvedCount\\\":0,\\\"solveKey\\\":\\\"default\\\"}}\",\"smallData\":\"{}\",\"solvable\":1,\"tags\":\"[\\\"VARIABLE\\\"]\",\"timestamp\":$(JSON.json(v.timestamp)),\"_version\":\"0.18.1\"}"

  v = Variable("x0", :Pose1, ["TEST_TAG"])
  @test JSON.json(v) == "{\"label\":\"x0\",\"dataEntry\":\"{}\",\"nstime\":\"0\",\"variableType\":\"IncrementalInference.ContinuousScalar\",\"dataEntryType\":\"{}\",\"ppeDict\":\"{}\",\"solverDataDict\":\"{\\\"default\\\":{\\\"vecval\\\":[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],\\\"dimval\\\":1,\\\"vecbw\\\":[0.0],\\\"dimbw\\\":1,\\\"BayesNetOutVertIDs\\\":[],\\\"dimIDs\\\":[0],\\\"dims\\\":1,\\\"eliminated\\\":false,\\\"BayesNetVertID\\\":\\\"_null\\\",\\\"separator\\\":[],\\\"variableType\\\":\\\"IncrementalInference.ContinuousScalar\\\",\\\"initialized\\\":false,\\\"infoPerCoord\\\":[0.0],\\\"ismargin\\\":false,\\\"dontmargin\\\":false,\\\"solveInProgress\\\":0,\\\"solvedCount\\\":0,\\\"solveKey\\\":\\\"default\\\"}}\",\"smallData\":\"{}\",\"solvable\":1,\"tags\":\"[\\\"TEST_TAG\\\"]\",\"timestamp\":$(JSON.json(v.timestamp)),\"_version\":\"0.18.1\"}"

  v = Variable("x0", :Pose2)
  @test JSON.json(v) == "{\"label\":\"x0\",\"dataEntry\":\"{}\",\"nstime\":\"0\",\"variableType\":\"RoME.Pose2\",\"dataEntryType\":\"{}\",\"ppeDict\":\"{}\",\"solverDataDict\":\"{\\\"default\\\":{\\\"vecval\\\":[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],\\\"dimval\\\":3,\\\"vecbw\\\":[0.0,0.0,0.0],\\\"dimbw\\\":3,\\\"BayesNetOutVertIDs\\\":[],\\\"dimIDs\\\":[0,1,2],\\\"dims\\\":3,\\\"eliminated\\\":false,\\\"BayesNetVertID\\\":\\\"_null\\\",\\\"separator\\\":[],\\\"variableType\\\":\\\"RoME.Pose2\\\",\\\"initialized\\\":false,\\\"infoPerCoord\\\":[0.0,0.0,0.0],\\\"ismargin\\\":false,\\\"dontmargin\\\":false,\\\"solveInProgress\\\":0,\\\"solvedCount\\\":0,\\\"solveKey\\\":\\\"default\\\"}}\",\"smallData\":\"{}\",\"solvable\":1,\"tags\":\"[\\\"VARIABLE\\\"]\",\"timestamp\":$(JSON.json(v.timestamp)),\"_version\":\"0.18.1\"}"

  v = Variable("x0", :Point2)
  @test JSON.json(v) == "{\"label\":\"x0\",\"dataEntry\":\"{}\",\"nstime\":\"0\",\"variableType\":\"RoME.Point2\",\"dataEntryType\":\"{}\",\"ppeDict\":\"{}\",\"solverDataDict\":\"{\\\"default\\\":{\\\"vecval\\\":[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],\\\"dimval\\\":2,\\\"vecbw\\\":[0.0,0.0],\\\"dimbw\\\":2,\\\"BayesNetOutVertIDs\\\":[],\\\"dimIDs\\\":[0,1],\\\"dims\\\":2,\\\"eliminated\\\":false,\\\"BayesNetVertID\\\":\\\"_null\\\",\\\"separator\\\":[],\\\"variableType\\\":\\\"RoME.Point2\\\",\\\"initialized\\\":false,\\\"infoPerCoord\\\":[0.0,0.0],\\\"ismargin\\\":false,\\\"dontmargin\\\":false,\\\"solveInProgress\\\":0,\\\"solvedCount\\\":0,\\\"solveKey\\\":\\\"default\\\"}}\",\"smallData\":\"{}\",\"solvable\":1,\"tags\":\"[\\\"VARIABLE\\\"]\",\"timestamp\":$(JSON.json(v.timestamp)),\"_version\":\"0.18.1\"}"

  v = Variable("x0", "RoME.Pose2")
  @test JSON.json(v) == "{\"label\":\"x0\",\"dataEntry\":\"{}\",\"nstime\":\"0\",\"variableType\":\"RoME.Pose2\",\"dataEntryType\":\"{}\",\"ppeDict\":\"{}\",\"solverDataDict\":\"{\\\"default\\\":{\\\"vecval\\\":[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],\\\"dimval\\\":3,\\\"vecbw\\\":[0.0,0.0,0.0],\\\"dimbw\\\":3,\\\"BayesNetOutVertIDs\\\":[],\\\"dimIDs\\\":[0,1,2],\\\"dims\\\":3,\\\"eliminated\\\":false,\\\"BayesNetVertID\\\":\\\"_null\\\",\\\"separator\\\":[],\\\"variableType\\\":\\\"RoME.Pose2\\\",\\\"initialized\\\":false,\\\"infoPerCoord\\\":[0.0,0.0,0.0],\\\"ismargin\\\":false,\\\"dontmargin\\\":false,\\\"solveInProgress\\\":0,\\\"solvedCount\\\":0,\\\"solveKey\\\":\\\"default\\\"}}\",\"smallData\":\"{}\",\"solvable\":1,\"tags\":\"[\\\"VARIABLE\\\"]\",\"timestamp\":$(JSON.json(v.timestamp)),\"_version\":\"0.18.1\"}"

end


text:
using FixedPointAcceleration
# Testing Error Evaluating Function
simple_vector_function(x) = [0.5*sqrt(x[1] + x[2]), 1.5*x[1] + 0.5*x[2]]
Inputs = [0.3,900]
fp = fixed_point(simple_vector_function, Inputs; Algorithm = Anderson)
# Inspecting this fp reveals an error after the 3rd iteration because
# Anderson tries to use a negative value for both x entries which results in
# the square root of a negative number. We can switch to simple
# iterations for a while to fix this.
fp.TerminationCondition_ == InvalidInputOrOutputOfIteration
fp.FailedEvaluation_.Error_ == ErrorExecutingFunction
fp = fixed_point(simple_vector_function, fp; Algorithm = Simple, MaxIter = Integer(7))
fp.TerminationCondition_ == ReachedMaxIter
fp = fixed_point(simple_vector_function, fp; Algorithm = Anderson)
fp.TerminationCondition_ == ReachedConvergenceThreshold

# Testing Input of NaN
fp = fixed_point(simple_vector_function, [NaN,900])
fp.FailedEvaluation_.Error_ == InputNAsDetected
# Testing Input of Inf
fp = fixed_point(simple_vector_function, [-Inf,900])
fp.FailedEvaluation_.Error_ == InputInfsDetected

# Testing Output of Nan
function funcfunc(x::Array{Float64,1})
    if abs(x[1] - 4.0) < 1e-12
        return Array{Float64,1}([NaN,4.0])
    end
    return sqrt.(x)
end
Inputs = [4.0,1.0]
fp = fixed_point(funcfunc, Inputs; Algorithm = Anderson)
fp.FailedEvaluation_.Error_ == OutputNAsDetected
# Testing Output of Missing
function funcfunc(x::Array{Float64,1})
    if abs(x[1] - 4.0) < 1e-12
        return [missing,4.0]
    end
    return sqrt.(x)
end
Inputs = [4.0,1.0]
fp = fixed_point(funcfunc, Inputs; Algorithm = Anderson)
fp.FailedEvaluation_.Error_ == OutputMissingsDetected
# Testing Output of Inf
function funcfunc(x::Array{Float64,1})
    if abs(x[1] - 4.0) < 1e-12
        return Array{Float64,1}([Inf,4.0])
    end
    return sqrt.(x)
end
Inputs = [4.0,1.0]
fp = fixed_point(funcfunc, Inputs; Algorithm = Anderson)
fp.FailedEvaluation_.Error_ == OutputInfsDetected
# Testing Output of wrong size
function funcfunc(x::Array{Float64,1})
    if abs(x[1] - 4.0) < 1e-12
        return Array{Float64,1}([5.0,4.0, 4.0])
    end
    return sqrt.(x)
end
Inputs = [4.0,1.0]
fp = fixed_point(funcfunc, Inputs; Algorithm = Anderson)
fp.FailedEvaluation_.Error_ == LengthOfOutputNotSameAsInput

# Testing Output of wrong type
function funcfunc(x::Array{Float64,1})
    return Array{Int,1}([5,4])
end
Inputs = [4.0,1.0]
fp = fixed_point(funcfunc, Inputs; Algorithm = Anderson)
fp.FailedEvaluation_.Error_ == FunctionIsNotTypeStable


text:
module MDPsTests

using MDPs

using Compat
using Base.Test
using FactCheck

include("fixtures.jl")

include("bellman.jl")
include("mdp.jl")
include("qfunction.jl")
include("random.jl")
include("reward.jl")
include("sparse_random.jl")
include("transition.jl")
include("utils.jl")

FactCheck.exitstatus()

end # module


text:
_berror() = throw(ArgumentError("Number of output images must be equal to the number of input images"))

inputvector(inputs::AbstractArray, args...) = obsview(inputs, args...)
inputvector(inputs::Tuple{Vararg{AbstractArray}}, args...) = obsview(inputs, args...)
@inline inputvector(inputs::AbstractVector{<:AbstractArray}, args...) = inputs
@inline inputvector(inputs::AbstractVector{<:Tuple{Vararg{AbstractArray}}}, args...) = inputs

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

"""
    augmentbatch!([resource], outputs, inputs, pipeline, [obsdim]) -> outputs

Apply the operations of the given `pipeline` to the images in
`inputs` and write the resulting images into `outputs`.

Both `outputs` and `inputs` have to contain the same number of images.
Each of these two variables can either be in the form of a higher
dimensional array, in the form of a vector of arrays for which
each vector element denotes an image.

```julia
# create five example observations of size 3x3
inputs = rand(3,3,5)
# create output arrays of appropriate shape
outputs = similar(inputs)
# transform the batch of images
augmentbatch!(outputs, inputs, FlipX() |> FlipY())
```

If one (or both) of the two parameters `outputs` and `inputs` is a
higher dimensional array, then the optional parameter `obsdim`
can be used specify which dimension denotes the observations
(defaults to `ObsDim.Last()`),

```julia
# create five example observations of size 3x3
inputs = rand(5,3,3)
# create output arrays of appropriate shape
outputs = similar(inputs)
# transform the batch of images
augmentbatch!(outputs, inputs, FlipX() |> FlipY(), ObsDim.First())
```

Similar to [`augment!`](@ref), it is also allowed for `outputs` and
`inputs` to both be tuples of the same length. If that is the case,
then each tuple element can be in any of the forms listed above.
This is useful for tasks such as image segmentation, where each
observations is made up of more than one image.

```julia
# create five example observations where each observation is
# made up of two conceptually linked 3x3 arrays
inputs = (rand(3,3,5), rand(3,3,5))
# create output arrays of appropriate shape
outputs = similar.(inputs)
# transform the batch of images
augmentbatch!(outputs, inputs, FlipX() |> FlipY())
```

The parameter `pipeline` can be a `AugmentorBase.Pipeline`, a tuple
of `AugmentorBase.Operation`, or a single `AugmentorBase.Operation`.

```julia
augmentbatch!(outputs, inputs, FlipX() |> FlipY())
augmentbatch!(outputs, inputs, (FlipX(), FlipY()))
augmentbatch!(outputs, inputs, FlipX())
```

The optional first parameter `resource` can either be `CPU1()`
(default) or `CPUThreads()`. In the later case the images will be
augmented in parallel. For this to make sense make sure that the
environment variable `JULIA_NUM_THREADS` is set to a reasonable
number so that `Threads.nthreads()` is greater than 1.

```julia
# transform the batch of images in parallel using multithreading
augmentbatch!(CPUThreads(), outputs, inputs, FlipX() |> FlipY())
```
"""
function augmentbatch!(
        outputs::Union{Tuple, AbstractArray},
        inputs::Union{Tuple, AbstractArray},
        pipeline,
        args...)
    augmentbatch!(CPU1(), outputs, inputs, pipeline, args...)
end

function augmentbatch!(
        r::AbstractResource,
        outputs::Union{Tuple, AbstractArray},
        inputs::Union{Tuple, AbstractArray},
        pipeline,
        obsdim = MLDataPattern.default_obsdim(outputs))
    augmentbatch!(r, inputvector(outputs, obsdim), inputvector(inputs, obsdim), pipeline)
    outputs
end

function augmentbatch!(
        ::CPU1,
        outputs::AbstractVector{<:Union{Tuple, AbstractArray}},
        inputs::AbstractVector{<:Union{Tuple, AbstractArray}},
        pipeline)
    length(outputs) == length(inputs) || _berror()
    for i in 1:length(outputs)
        augment!(outputs[i], inputs[i], pipeline)
    end
    outputs
end

function augmentbatch!(
        ::CPUThreads,
        outputs::AbstractVector{<:Union{Tuple, AbstractArray}},
        inputs::AbstractVector{<:Union{Tuple, AbstractArray}},
        pipeline)
    length(outputs) == length(inputs) || _berror()
    Threads.@threads for i in 1:length(outputs)
        augment!(outputs[i], inputs[i], pipeline)
    end
    outputs
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 4155e7d0-a9f6-11eb-1cf4-ff8e390713db
using SpectralVis

# ╔═╡ 474e20d1-fce7-4d63-95b5-511af903b2a9
using Plots

# ╔═╡ 8af080d9-c436-4de7-8443-71b8b658e7dd
function spectralmix(c1, c2)
	reflectance_spec(c1.λ,(c1.r.+c2.r)/2)
end

# ╔═╡ 0a4f8542-87a6-4e66-a895-94eaded147c0
set_specenv()

# ╔═╡ faeb5ddd-37ed-47b3-a2bf-062a2635fa98
SPECENV

# ╔═╡ 87eb6bed-da99-46cc-a1fb-30a561026331
begin
	blue=macbeth(:Macbeth_chart2,24)
	red=macbeth(:Macbeth_chart2,19)
	mix=spectralmix(red,blue)
	cmatch = SPECENV.cmf
	rspl = spline3(red.λ,red.r,:natural)
	bspl = spline3(blue.λ,blue.r,:natural)
	pspl = spline3(mix.λ,mix.r,:natural)
	re = extrap(rspl,390:830,:linear)
	be = extrap(bspl,390:830,:linear)
	mixe = extrap(pspl,390:830,:linear)
	ri = interp(re,390:1:830)
	bi = interp(be,390:1:830)
	mixi = interp(mixe,390:1:830)
	dill=D_series_illuminant(6500)
	dswp=D_series_whitepoint(6500)
	redcol=Colors.convert(Lab,dill*RSpec(ri[1],ri[2])*SPECENV.cmf,dswp)
	blucol=Colors.convert(Lab,dill*RSpec(bi[1],bi[2])*SPECENV.cmf,dswp)
	mixcol=Colors.convert(Lab,dill*RSpec(mixi[1],mixi[2])*SPECENV.cmf,dswp)

	scatter([blucol.a mixcol.a redcol.a],[blucol.b mixcol.b redcol.b], color=[blucol mixcol redcol],markersize=10,legend=:outertopright,xlimits=(-55,55),ylimits=(-55,55),aspect_ratio=1)
end

# ╔═╡ bf708eec-143a-4b6a-9f6d-f163134e0724
blucol.l

# ╔═╡ cd316ddd-3e70-4769-8686-b0fd7f034ab4
dill*RSpec(ri[1],ri[2])*SPECENV.cmf

# ╔═╡ 7b21efd3-16ff-4a9c-a6ca-4c8902b0c7a1
dswp

# ╔═╡ deecc0c1-cc04-49d5-aa99-8708f77f6ca1
[blucol,mixcol,redcol]

# ╔═╡ Cell order:
# ╠═4155e7d0-a9f6-11eb-1cf4-ff8e390713db
# ╠═474e20d1-fce7-4d63-95b5-511af903b2a9
# ╠═8af080d9-c436-4de7-8443-71b8b658e7dd
# ╠═0a4f8542-87a6-4e66-a895-94eaded147c0
# ╠═faeb5ddd-37ed-47b3-a2bf-062a2635fa98
# ╠═87eb6bed-da99-46cc-a1fb-30a561026331
# ╠═bf708eec-143a-4b6a-9f6d-f163134e0724
# ╠═cd316ddd-3e70-4769-8686-b0fd7f034ab4
# ╠═7b21efd3-16ff-4a9c-a6ca-4c8902b0c7a1
# ╠═deecc0c1-cc04-49d5-aa99-8708f77f6ca1


text:
typealias CudaDims Union(Integer, (Integer,), (Integer, Integer), (Integer, Integer, Integer))

const CudaExtent = rt.cudaExtent
CudaExtent(dims::CudaDims) = CudaExtent(get_size_x(dims), get_size_y(dims), get_size_z(dims))
CudaExtent{T}(::Type{T}, dims::CudaDims) = CudaExtent(get_size_x(dims)*sizeof(T), get_size_y(dims), get_size_z(dims))

const CudaPos = rt.cudaPos
CudaPos(dims::CudaDims) = CudaPos(get_size_x(dims)-1, get_size_y(dims)-1, get_size_z(dims)-1)
CudaPos{T}(::Type{T}, dims::CudaDims) = CudaPos((get_size_x(dims)-1)*sizeof(T), get_size_y(dims)-1, get_size_z(dims)-1)

get_size_x(g::Integer) = g
get_size_x(g::(Integer,)) = g[1]
get_size_x(g::(Integer, Integer)) = g[1]
get_size_x(g::(Integer, Integer, Integer)) = g[1]

get_size_y(g::Integer) = 1
get_size_y(g::(Integer,)) = 1
get_size_y(g::(Integer, Integer)) = g[2]
get_size_y(g::(Integer, Integer, Integer)) = g[2]

get_size_z(g::Integer) = 1
get_size_z(g::(Integer,)) = 1
get_size_z(g::(Integer, Integer)) = 1
get_size_z(g::(Integer, Integer, Integer)) = g[3]

# dim3 and utility functions
const Dim3 = rt.dim3

Dim3(dim::CudaDims) = Dim3(get_dim_x(dim), get_dim_y(dim), get_dim_z(dim))
Dim3(dims::Integer...) = Dim3(dims)

get_dim_x(g::CudaDims) = safe_uint32(get_size_x(g))
get_dim_y(g::CudaDims) = safe_uint32(get_size_y(g))
get_dim_z(g::CudaDims) = safe_uint32(get_size_z(g))

function safe_uint32(i::Integer)
    ui = uint32(i)
    if ui < i
        error("Invalid integer size")
    end
    ui
end


text:
function run_investor_iteration(investor::Investor,
                                 active_projects::Vector{Project},
                                 iteration_year::Int64,
                                 yearly_horizon::Int64,
                                 simulation_years::Int64,
                                 start_year::Int64,
                                 capacity_forward_years::Int64,
                                 sys_data_dir::String,
                                 sys_UC::Union{Nothing, PSY.System},
                                 services::Vector{PSY.Service},
                                 solver::JuMP.MOI.OptimizerWithAttributes
                            )

    investor_dir = get_data_dir(investor)
    projects = get_projects(investor)

    active_projects_copy = deepcopy(active_projects)

    market_names = get_markets(investor)

    # Create empty market prices struct
    market_prices = MarketPrices()

    option_projects = get_options(investor)
    max_new_options = Dict(get_name(project) => 0 for project in option_projects)

    scenarios = get_scenario_data(get_forecast(investor))

    for scenario in scenarios
        scenario_name = get_name(scenario)
        output_file = joinpath(investor_dir, "expected_market_data", "$(scenario_name)_year_$(iteration_year).jld2")
        expected_data = FileIO.load(output_file)

        set_energy_price!(market_prices, scenario_name, expected_data["energy_price"])

        if in(:ReserveUp, market_names)
            set_reserveup_price!(market_prices, scenario_name, expected_data["reserve_up_price"])
        end

        if in(:ReserveDown, market_names)
            set_reservedown_price!(market_prices, scenario_name, expected_data["reserve_down_price"])
        end

        if in(:Capacity, market_names)
            set_capacity_price!(market_prices, scenario_name, expected_data["capacity_price"])
        end

        if in(:REC, market_names)
            set_rec_price!(market_prices, scenario_name, expected_data["rec_price"])
        end

        for project in projects
            if in(get_name(project), get_name.(active_projects_copy))
                update_capacity_factors!(project, scenario_name, expected_data["capacity_factors"])
            end
            if in(get_name(project), get_name.(active_projects_copy))
                update_capacity_accepted_perc!(project, scenario_name, expected_data["capacity_accepted_perc"])
            end

        end

         max_new_options = update_max_new_options!(max_new_options, expected_data["new_options"], option_projects)
    end

    set_market_prices!(investor, market_prices)

    retire_unprofitable!(investor,
                         sys_UC,
                         sys_data_dir,
                         iteration_year,
                         yearly_horizon,
                         simulation_years,
                         capacity_forward_years,
                         solver)

    make_investments!(investor,
                      max_new_options,
                      iteration_year,
                      yearly_horizon,
                      simulation_years,
                      capacity_forward_years,
                      solver)

    for (i, project) in enumerate(projects)
        start_construction!(projects,
                            i,
                            project,
                            iteration_year)

        finish_construction!(projects,
                            i,
                            project,
                            sys_UC,
                            services,
                            sys_data_dir,
                            iteration_year,
                            start_year)

        update_lifecycle!(project,
                          iteration_year,
                          simulation_years)
    end

    return
end


text:
export iauS06a
"""
The CIO locator s, positioning the Celestial Intermediate Origin on
the equator of the Celestial Intermediate Pole, using the IAU 2006
precession and IAU 2000A nutation models.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   date1,date2  double    TT as a 2-part Julian Date (Note 1)

Returned (function value):
               double    the CIO locator s in radians (Note 2)

Notes:

   1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   2. The CIO locator s is the difference between the right ascensions
   of the same point in two systems.  The two systems are the GCRS
   and the CIP,CIO, and the point is the ascending node of the
   CIP equator.  The CIO locator s remains a small fraction of
   1 arcsecond throughout 1900-2100.

   3. The series used to compute s is in fact for s+XY/2, where X and Y
   are the x and y components of the CIP unit vector;  this series is
   more compact than a direct series for s would be.  The present
   function uses the full IAU 2000A nutation model when predicting
   the CIP position.

Called:
   iauPnm06a    classical NPB matrix, IAU 2006/2000A
   iauBpn2xy    extract CIP X,Y coordinates from NPB matrix
   iauS06       the CIO locator s, given X,Y, IAU 2006

References:

   Capitaine, N., Chapront, J., Lambert, S. and Wallace, P.,
   "Expressions for the Celestial Intermediate Pole and Celestial
   Ephemeris Origin consistent with the IAU 2000A precession-
   nutation model", Astron.Astrophys. 400, 1145-1154 (2003)

   n.b. The celestial ephemeris origin (CEO) was renamed "celestial
         intermediate origin" (CIO) by IAU 2006 Resolution 2.

   Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855

   McCarthy, D. D., Petit, G. (eds.), 2004, IERS Conventions (2003),
   IERS Technical Note No. 32, BKG

   Wallace, P.T. & Capitaine, N., 2006, Astron.Astrophys. 459, 981

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# double iauS06a(double date1, double date2)
function iauS06a(date1::Real, date2::Real)
   return ccall((:iauS06a, libsofa_c), Cdouble,
                   (Cdouble, Cdouble),
                   convert(Float64, date1), convert(Float64, date2))
end


text:
function linsolve(operator, b, x₀, alg::GMRES, a₀::Number = 0, a₁::Number = 1)
    # Initial function operation and division defines number type
    y₀ = apply(operator, x₀)
    T = typeof(dot(b, y₀)/norm(b)*one(a₀)*one(a₁))
    α₀ = convert(T, a₀)::T
    α₁ = convert(T, a₁)::T
    # Continue computing r = b - a₀ * x₀ - a₁ * operator(x₀)
    r = copyto!(similar(b, T), b)
    r = α₀ == 0 ? r : axpy!(-α₀, x₀, r)
    r = axpy!(-α₁, y₀, r)
    x = copyto!(similar(r), x₀)
    β = norm(r)
    S = typeof(β)

    # Algorithm parameters
    maxiter = alg.maxiter
    krylovdim = alg.krylovdim
    tol = convert(S, alg.rtol == 0 ? alg.atol : max(alg.atol, alg.rtol*norm(b)))::S

    # Check for early return
    β < tol && return (x, ConvergenceInfo(1, r, β, 0, 1))

    # Initialize data structures
    y = Vector{T}(undef, krylovdim+1)
    gs = Vector{Givens{T}}(undef, krylovdim)
    R = fill(zero(T), (krylovdim,krylovdim))
    numiter = 0
    numops = 1 # operator has been applied once to determine T

    iter = ArnoldiIterator(operator, r, alg.orth)
    fact = initialize(iter)
    numops += 1 # start applies operator once
    while numiter < maxiter # restart loop
        numiter += 1
        y[1] = β
        k = 1
        H = rayleighquotient(fact)
        R[1,1] = α₀ + α₁ * H[1,1]
        gs[1], R[1,1] = givens(R[1,1], α₁*normres(fact), 1, 2)
        y[2] = zero(T)
        lmul!(gs[1], y)
        β = convert(S, abs(y[2]))
        # info("iter $numiter, step $k : normres = $β")

        while β > tol && length(fact) < krylovdim # inner arnoldi loop
            fact = expand!(iter, fact)
            numops += 1 # expand! applies the operator once
            k = length(fact)
            H = rayleighquotient(fact)

            # copy Arnoldi Hessenberg matrix into R
            @inbounds begin
                for i=1:k-1
                    R[i,k] = α₁ * H[i,k]
                end
                R[k,k] = α₀ + α₁ * H[k,k]
            end

            # Apply Givens rotations
            Rk = view(R, :, k)
            @inbounds for i=1:k-1
                lmul!(gs[i], Rk)
            end
            gs[k], R[k,k] = givens(R[k,k], α₁*normres(fact), k, k+1)

            # Apply Givens rotations to right hand side
            y[k+1] = zero(T)
            lmul!(gs[k], y)

            # New error
            β = convert(S, abs(y[k+1]))
            # info("iter $numiter, step $k : normres = $β")
        end
        # info("iter $numiter, finished at step $k : normres = $β")

        # Solve upper triangular system
        ldiv!(UpperTriangular(R), y, 1:k)

        # Update x
        V = basis(fact)
        @inbounds for i = 1:k
            x = axpy!(y[i], V[i], x)
        end

        if β > tol
            # Recompute residual without reevaluating operator
            w = residual(fact)
            push!(V, rmul!(w, 1/normres(fact)))
            for i = 1:k
                rmul!(V, gs[i]')
            end
            r = mul!(r, y[k+1], V[k+1])
        else
            # Recompute residual and its norm explicitly, to ensure that no
            # numerical errors have accumulated
            r = copyto!(r, b)                      # r = b
            r = α₀ == 0 ? r : axpy!(-α₀, x, r)     #      - α₀ * x
            r = axpy!(-α₁, apply(operator, x), r)  #      - α₁ * operator(x)
            numops += 1
            β = norm(r)
            β < tol && return (x, ConvergenceInfo(1, r, β, numiter, numops))
        end

        # Restart Arnoldi factorization with new r
        iter = ArnoldiIterator(operator, r, alg.orth)
        fact = initialize!(iter, fact)
    end
    return (x, ConvergenceInfo(0, r, β, numiter, numops))
end


text:
## This file contains libblastrampoline-specific APIs

# Keep these in sync with `src/libblastrampoline_internal.h`
struct lbt_library_info_t
    libname::Cstring
    handle::Ptr{Cvoid}
    suffix::Cstring
    active_forwards::Ptr{UInt8}
    interface::Int32
    f2c::Int32
end
const LBT_INTERFACE_LP64    = 32
const LBT_INTERFACE_ILP64   = 64
const LBT_INTERFACE_UNKNOWN = -1
const LBT_INTERFACE_MAP = Dict(
    LBT_INTERFACE_LP64    => :lp64,
    LBT_INTERFACE_ILP64   => :ilp64,
    LBT_INTERFACE_UNKNOWN => :unknown,
)

const LBT_F2C_PLAIN         =  0
const LBT_F2C_REQUIRED      =  1
const LBT_F2C_UNKNOWN       = -1
const LBT_F2C_MAP = Dict(
    LBT_F2C_PLAIN    => :plain,
    LBT_F2C_REQUIRED => :required,
    LBT_F2C_UNKNOWN  => :unknown,
)

struct LBTLibraryInfo
    libname::String
    handle::Ptr{Cvoid}
    suffix::String
    active_forwards::Vector{UInt8}
    interface::Symbol
    f2c::Symbol

    function LBTLibraryInfo(lib_info::lbt_library_info_t, num_exported_symbols::UInt32)
        return new(
            unsafe_string(lib_info.libname),
            lib_info.handle,
            unsafe_string(lib_info.suffix),
            unsafe_wrap(Vector{UInt8}, lib_info.active_forwards, div(num_exported_symbols,8)+1),
            LBT_INTERFACE_MAP[lib_info.interface],
            LBT_F2C_MAP[lib_info.f2c],
        )
    end
end

struct lbt_config_t
    loaded_libs::Ptr{Ptr{lbt_library_info_t}}
    build_flags::UInt32
    exported_symbols::Ptr{Cstring}
    num_exported_symbols::UInt32
end
const LBT_BUILDFLAGS_DEEPBINDLESS = 0x01
const LBT_BUILDFLAGS_F2C_CAPABLE  = 0x02
const LBT_BUILDFLAGS_MAP = Dict(
    LBT_BUILDFLAGS_DEEPBINDLESS => :deepbindless,
    LBT_BUILDFLAGS_F2C_CAPABLE => :f2c_capable,
)

struct LBTConfig
    loaded_libs::Vector{LBTLibraryInfo}
    build_flags::Vector{Symbol}
    exported_symbols::Vector{String}

    function LBTConfig(config::lbt_config_t)
        # Decode OR'ed flags into a list of names
        build_flag_names = Symbol[]
        for (flag, name) in LBT_BUILDFLAGS_MAP
            if config.build_flags & flag != 0x00
                push!(build_flag_names, name)
            end
        end

        # Load all exported symbol names
        exported_symbols = String[]
        for sym_idx in 1:config.num_exported_symbols
            str_ptr = unsafe_load(config.exported_symbols, sym_idx)
            if str_ptr != C_NULL
                push!(exported_symbols, unsafe_string(str_ptr))
            else
                @error("NULL string in lbt_config.exported_symbols[$(sym_idx)]")
            end
        end

        # Unpack library info structures
        libs = LBTLibraryInfo[]
        idx = 1
        lib_ptr = unsafe_load(config.loaded_libs, idx)
        while lib_ptr != C_NULL
            push!(libs, LBTLibraryInfo(unsafe_load(lib_ptr), config.num_exported_symbols))

            idx += 1
            lib_ptr = unsafe_load(config.loaded_libs, idx)
        end
        return new(
            libs,
            build_flag_names,
            exported_symbols,
        )
    end
end

function lbt_get_config()
    config_ptr = ccall((:lbt_get_config, libblastrampoline), Ptr{lbt_config_t}, ())
    return LBTConfig(unsafe_load(config_ptr))
end

function lbt_get_num_threads()
    return ccall((:lbt_get_num_threads, libblastrampoline), Int32, ())
end

function lbt_set_num_threads(nthreads)
    return ccall((:lbt_set_num_threads, libblastrampoline), Cvoid, (Int32,), nthreads)
end

function lbt_forward(path; clear::Bool = false, verbose::Bool = false)
    ccall((:lbt_forward, libblastrampoline), Int32, (Cstring, Int32, Int32), path, clear ? 1 : 0, verbose ? 1 : 0)
end

function lbt_set_default_func(addr)
    return ccall((:lbt_set_default_func, libblastrampoline), Cvoid, (Ptr{Cvoid},), addr)
end

function lbt_get_default_func()
    return ccall((:lbt_get_default_func, libblastrampoline), Ptr{Cvoid}, ())
end

#=
Don't define footgun API (yet)

function lbt_get_forward(symbol_name, interface, f2c = LBT_F2C_PLAIN)
    return ccall((:lbt_get_forward, libblastrampoline), Ptr{Cvoid}, (Cstring, Int32, Int32), symbol_name, interface, f2c)
end

function lbt_set_forward(symbol_name, addr, interface, f2c = LBT_F2C_PLAIN; verbose::Bool = false)
    return ccall((:lbt_set_forward, libblastrampoline), Int32, (Cstring, Ptr{Cvoid}, Int32, Int32, Int32), symbol_name, addr, interface, f2c, verbose ? 1 : 0)
end
=#


text:
#!/usr/bin/env julia

import Celeste.DeterministicVI: elbo, ElboArgs
using DiffBase

include(string(Pkg.dir("Celeste"), "/test/Synthetic.jl"))
include(string(Pkg.dir("Celeste"), "/test/SampleData.jl"))

function main()
    srand(1)
    println("Loading data.")

    S = 100
    blob, ea, body = SampleData.gen_n_body_dataset(S)
    ea = ElboArgs(ea.images, ea.vp, ea.patches, ea.active_sources;
                  psf_K=ea.psf_K,
                  calculate_gradient=true,
                  calculate_hessian=false)

    param_length = length(Celeste.Model.CanonicalParams)
    kl_source = Celeste.SensitiveFloats.SensitiveFloat{Float64}(param_length, 1, true, false)

    println("Warm-up / compiling.")
    # do a trial run first, so we don't profile/time compling the code
    elbo(ea, kl_source)
    Profile.clear_malloc_data()

    println("Calculating ELBO and gradient.")
    if isempty(ARGS)
        # let's time it without any overhead from profiling
        @time elbo(ea, kl_source)
    elseif ARGS[1] == "--profile"
        Profile.init(delay=0.001)
        @profile elbo(ea, kl_source)
        Profile.print(format=:flat, sortedby=:count)
    end
end

main()


text:
#
# SubGradient Plan
#
export SubGradientMethodOptions, SubGradientProblem
export getCost, getSubGradient
export getStepsize
#
# Problem
@doc doc"""
    SubGradientProblem <: Problem
A structure to store information about a subgradient based optimization problem

# Fields
* `M` – a [`Manifold`](@ref)
* `costFunction` – the function $F$ to be minimized
* `subGradient` – a function returning a subgradient $\partial F$ of $F$
"""
mutable struct SubGradientProblem{mT <: Manifold} <: Problem
    M::mT
    costFunction::Function
    subGradient::Function
end
"""
    getSubGradient(p,x)

evaluate the gradient of a [`SubGradientProblem`](@ref)` p` at the [`MPoint`](@ref) `x`.
"""
getSubGradient(p::P,x::MP) where {P <: SubGradientProblem{M} where M <: Manifold, MP <: MPoint} = p.subGradient(x)
"""
    SubGradientMethodOptions <: Options
stories option values for a [`subGradientMethod`](@ref) solver

# Fields
* `retraction` – the retration to use within
* `stepsize` – a [`Stepsize`](@ref)
* `stop` – a [`StoppingCriterion`](@ref)
* `x` – (initial or current) value the algorithm is at
* `optimalX` – optimal value
"""
mutable struct SubGradientMethodOptions{P,T} <: Options where {P <: MPoint, T <: TVector}
    retraction::Function
    stepsize::Stepsize
    stop::StoppingCriterion
    x::P
    xOptimal::P
    ∂::T
    SubGradientMethodOptions{P,T}(x::P,sC::StoppingCriterion,s::Stepsize,retr::Function=exp) where {P <: MPoint, T <: TVector} = (
        o = new{P,T}(); o.x = x; o.xOptimal = x;
        o.stepsize = s; o.retraction = retr;
        o.stop = sC;
        return o
    )
end
SubGradientMethodOptions(x::P,sC::StoppingCriterion,s::Stepsize,retr::Function=exp) where {P <: MPoint} = SubGradientMethodOptions{P,typeofTVector(x)}(x,sC,s,retr)


text:
using FastGaussQuadrature

"""
    get_quadpoints(n)

Compute quadrature points and weights for Gaussian quadrature
of order `n`.
The points (and thus the corresponding weights) are normalized to the range
``[0.0, 1.0]``.

Return a tuple of `points, weights`.
"""
function get_quadpoints(n)
    x, w = gausslegendre(n)
    # They are scaled in [-1.0, 1.0]
    # and we need [0.0, 1.0]
    (x .+ 1)./2, w ./ 2
end


text:
@testset "grm" begin
    vcf = "test.08Jun17.d8b.vcf.gz"
    isfile(vcf) || download("http://faculty.washington.edu/browning/beagle/test.08Jun17.d8b.vcf.gz",
        abspath(joinpath(dirname(pathof(VCFTools)), "..", "test/$vcf")))
    # standard GRM method
    Φgrm = grm(vcf, method=:GRM)
    n = nsamples(vcf)
    p = nrecords(vcf)
    @test size(Φgrm) == (n, n)
    @test issymmetric(Φgrm)
    @test eigmin(Φgrm) > -1e-8
    # MoM estimates
    Φmom = grm(vcf, method=:MoM)
    @test size(Φmom) == (n, n)
    @test issymmetric(Φmom)
    # Robust GRM
    Φrbs = grm(vcf, method=:Robust)
    @test size(Φrbs) == (n, n)
    @test issymmetric(Φrbs)
    @test eigmin(Φrbs) > -1e-8
    # test scale_missing
    Random.seed!(2020)
    masked_vcf = "masked.test.08Jun17.d8b.vcf.gz"
    mask_gt(vcf, bitrand(p, n), des=masked_vcf) # create vcf file with missings
    Φrbs_sm = grm(masked_vcf, method=:Robust, scale_missing=true)
    @test issymmetric(Φrbs_sm)
    @test eigmin(Φrbs_sm) > -1e-8
    @test size(Φrbs_sm) == (n, n)
end


text:
@testset "FBD.used_bits" begin
    @test FBD.used_bits(5) == 3
    @test FBD.used_bits(0) == 0
end

@testset "FBD.kronecker_product" begin
    a = [0 2; 3 4]
    p0, p1, p2 = FBD.kronecker_power.([a], 0:2)
    @test p0 == reshape([1], (1,1))
    @test p1 === a
    @test p2 == [0 0 0 4; 0 0 6 8; 0 6 0 8; 9 12 12 16]
    @test reduce(FBD.kronecker_product, [p0, p1, p0, p0, p1]) == p2
end

max_stars_repo_name: LilithHafner/FBD

max_issues_repo_name: LilithHafner/FBD

max_forks_repo_name: LilithHafner/FBD


text:
using Documenter
using JSON
using ThreadsX

function should_push_preview(event_path = get(ENV, "GITHUB_EVENT_PATH", nothing))
    event_path === nothing && return false
    event = JSON.parsefile(event_path)
    pull_request = get(event, "pull_request", nothing)
    pull_request === nothing && return false
    labels = [x["name"] for x in pull_request["labels"]]
    # https://developer.github.com/v3/activity/events/types/#pullrequestevent
    yes = "push_preview" in labels
    if yes
        @info "Trying to push preview as label `push_preview` is specified." labels
    else
        @info "Not pushing preview as label `push_preview` is not specified." labels
    end
    return yes
end

makedocs(;
    modules=[ThreadsX],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/tkf/ThreadsX.jl/blob/{commit}{path}#L{line}",
    sitename="ThreadsX.jl",
    authors="Takafumi Arakaki <aka.tkf@gmail.com>",
)

deploydocs(;
    repo="github.com/tkf/ThreadsX.jl",
    push_preview = should_push_preview(),
)