text:
module SharedControl

using MathProgBase # for checkFeasibility!

include("CaseModule.jl")
using .CaseModule

export
      initializeSharedControl!,
      checkFeasibility!,
      sharedControl!,
      getPlantData!,
      sendOptData!,
      ExternalModel

type ExternalModel  # communication
  s1
  s2
  status   # to pass result of optimization to Matlab
  sa_sample
  runJulia # a Bool (in Int form) to indicate to run julia comming from Matlab
  numObs   # number of obstacles
  SA       # last two drivers steering angle
  SA_first # first order approximation of drivers commands
  UX       # vehicle speed
  X0       # vehicle states
  Xobs
  Yobs
  A
  B
  t0  # current time
  t_opt
  wsa
  feasible # a bool for the feasibility of the current driver commands
  infeasible_counter
  infeasible_counter_max
  AL
  AU
  dt  # this is the time step in simulink
end

function ExternalModel()
  ExternalModel(Any,
                Any,
                1.0,
                [],
                1,
                3,
                [0.0,0.0],
                [],
                10.0,
                [],
                [],
                [],
                [],
                [],
                0.0,
                0.0,
                1.0,
                false,
                100,
                10,
                [],
                [],
                0.01)
end

"""
mdl,n,r,params,x,d=initializeSharedControl!(c)
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 1/27/2017, Last Modified: 4/4/2017 \n
--------------------------------------------------------------------------------------\n
"""
function initializeSharedControl!(c)
    if c.m.model!=:ThreeDOFv2; error(" this is for ThreeDOFv2 ! ") end

    if c.m.activeSafety
      pa=Vpara(x_min=c.m.Xlims[1],x_max=c.m.Xlims[2],y_min=c.m.Ylims[1],y_max=c.m.Ylims[2],sr_min=-0.18,sr_max=0.18);
    else
      pa=Vpara(x_min=c.m.Xlims[1],x_max=c.m.Xlims[2],y_min=c.m.Ylims[1],y_max=c.m.Ylims[2],sr_min=-0.10,sr_max=0.10); #0.12
    end
    @unpack_Vpara pa;

    XF=[  NaN, NaN,   NaN, NaN,     NaN,    NaN,    NaN, NaN];
  #  XL=[x_min, y_min, NaN, NaN, psi_min, sa_min, c.m.UX, 0.0];
  #  XU=[x_max, y_max, NaN, NaN, psi_max, sa_max, c.m.UX, 0.0];
    XL=[NaN, NaN, NaN, NaN, psi_min, sa_min, c.m.UX, 0.0];
    XU=[NaN, NaN, NaN, NaN, psi_max, sa_max, c.m.UX, 0.0];
    CL = [sr_min, 0.0]; CU = [sr_max, 0.0];
    n=define(ThreeDOFv2;numStates=8,numControls=2,X0=copy(c.m.X0),XF=XF,XL=XL,XU=XU,CL=CL,CU=CU);

    # variable names
             # 1  2  3  4  5    6   7   8
    names = [:x,:y,:v,:r,:psi,:sa,:ux,:ax];
    descriptions = ["X (m)","Y (m)","Lateral Velocity (m/s)", "Yaw Rate (rad/s)","Yaw Angle (rad)", "Steering Angle (rad)", "Longitudinal Velocity (m/s)", "Longitudinal Acceleration (m/s^2)"];
    stateNames!(n,names,descriptions);
             # 1    2
    names = [:sr,:jx];
    descriptions = ["Steering Rate (rad/s)","Longitudinal Jerk (m/s^3)"];
    controlNames!(n,names,descriptions);

    # configure problem
    configure!(n,Nck=c.m.Nck;(:integrationScheme => :lgrExplicit),(:finalTimeDV => false),(:tf => c.m.tp))
    mdl=defineSolver!(n,c);

    # define tolerances
    XF_tol=[NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN];
    X0_tol=[0.05,0.05,0.005,0.05,0.01,0.003,NaN,NaN];  # TODO BE CAREFUL HERE!!
    defineTolerances!(n;X0_tol=X0_tol,XF_tol=XF_tol);

    # add parameters
    @NLparameter(mdl, ux_param==c.m.UX);                       # inital vehicle speed
    @NLparameter(mdl, sa_param[i=1:n.numStatePoints]==0.0);    # driver steering command
    @NLparameter(mdl, sr_param==0.0);                          # driver steering rate command

    veh_params=[ux_param,sa_param,sr_param];

    # obstacles
    #Q=size(c.o.A)[1]; # number of obstacles TODO update these based off of LiDAR data
    Q=1;
    if Q>3; error("for activeSafety mode the code was modified to deal with only three obstales when checkFeasibility is running \n");end
    @NLparameter(mdl, a[i=1:Q] == c.o.A[i]);
    @NLparameter(mdl, b[i=1:Q] == c.o.B[i]);
    @NLparameter(mdl, X_0[i=1:Q] == c.o.X0[i]);
    @NLparameter(mdl, Y_0[i=1:Q] == c.o.Y0[i]);
    obs_params=[a,b,X_0,Y_0];

    x=ExternalModel(); # setup data with Matlab

    # setup mpc parameters
    initializeMPC!(n;FixedTp=c.m.FixedTp,PredictX0=c.m.PredictX0,tp=c.m.tp,tex=copy(c.m.tex),max_iter=c.m.mpc_max_iter);
    x.X0=copy(c.m.X0);

    # define ocp
    s=Settings(;save=false,MPC=true,evalConstraints=false);
    r=OCPdef!(mdl,n,s,[pa,ux_param]);  # need pa out of params -> also need speed for c.m.model==:ThreeDOFv1

    # define objective function
    obj=0;
    if c.m.followPath
      # follow the path -> min((X_path(Yt)-Xt)^2)
      if c.t.func==:poly
        path_obj=@NLexpression(mdl,c.w.path*sum(  (  (c.t.a[1] + c.t.a[2]*r.x[(i+1),2] + c.t.a[3]*r.x[(i+1),2]^2 + c.t.a[4]*r.x[(i+1),2]^3 + c.t.a[5]*r.x[(i+1),2]^4) - r.x[(i+1),1]  )^2 for i in 1:n.numStatePoints-1)  );
      elseif c.t.func==:fourier
        path_obj=@NLexpression(mdl,c.w.path*sum(  (  (c.t.a[1]*sin(c.t.b[1]*r.x[(i+1),1]+c.t.c[1]) + c.t.a[2]*sin(c.t.b[2]*r.x[(i+1),1]+c.t.c[2]) + c.t.a[3]*sin(c.t.b[3]*r.x[(i+1),1]+c.t.c[3]) + c.t.a[4]*sin(c.t.b[4]*r.x[(i+1),1]+c.t.c[4]) + c.t.a[5]*sin(c.t.b[5]*r.x[(i+1),1]+c.t.c[5]) + c.t.a[6]*sin(c.t.b[6]*r.x[(i+1),1]+c.t.c[6]) + c.t.a[7]*sin(c.t.b[7]*r.x[(i+1),1]+c.t.c[7]) + c.t.a[8]*sin(c.t.b[8]*r.x[(i+1),1]+c.t.c[8])+c.t.y0) - r.x[(i+1),2]  )^2 for i in 1:n.numStatePoints-1)  );
      end
      obj=path_obj;
    end

    if c.m.followDriver
       driver_obj=integrate!(mdl,n,r.x[:,6];C=c.w.driver,D=copy(sa_param),(:variable=>:control),(:integrand=>:squared),(:integrandAlgebra=>:subtract));
       obj=@NLexpression(mdl,obj+driver_obj);
    end
    sr_obj=integrate!(mdl,n,r.u[:,1];C=c.w.sr,(:variable=>:control),(:integrand=>:squared));     # minimize steering rate

    @NLobjective(mdl, Min, obj+sr_obj);

    # constraint steering rate
    sr_tol=0.005; #0.01
    sr_conL=@NLconstraint(mdl, [i=1],r.u[i,1] <=  (sr_param+sr_tol));
    newConstraint!(r,sr_conL,:sr_conL);
    sr_conU=@NLconstraint(mdl, [i=1],-r.u[i,1] <= -(sr_param-sr_tol));
    newConstraint!(r,sr_conU,:sr_conU);

    #sr_con=@NLconstraint(mdl,[i=1],r.u[i,1]==sr_param);
    #newConstraint!(r,sr_con,:sr_con);

    # obstacle postion after the intial postion
    @NLexpression(mdl, X_obs[j=1:Q,i=1:n.numStatePoints], X_0[j])
    @NLexpression(mdl, Y_obs[j=1:Q,i=1:n.numStatePoints], Y_0[j])

    # constraint position
    j=1;obs_con1=@NLconstraint(mdl, [i=3:n.numStatePoints-1], 1 <= ((r.x[(i+1),1]-X_obs[j,i])^2)/((a[j]+c.m.sm)^2) + ((r.x[(i+1),2]-Y_obs[j,i])^2)/((b[j]+c.m.sm)^2));
    newConstraint!(r,obs_con1,:obs_con1);
    #j=2;obs_con2=@NLconstraint(mdl, [i=1:n.numStatePoints-1], 1 <= ((r.x[(i+1),1]-X_obs[j,i])^2)/((a[j]+c.m.sm)^2) + ((r.x[(i+1),2]-Y_obs[j,i])^2)/((b[j]+c.m.sm)^2));
    #newConstraint!(r,obs_con2,:obs_con2);
    #j=3;obs_con3=@NLconstraint(mdl, [i=1:n.numStatePoints-1], 1 <= ((r.x[(i+1),1]-X_obs[j,i])^2)/((a[j]+c.m.sm)^2) + ((r.x[(i+1),2]-Y_obs[j,i])^2)/((b[j]+c.m.sm)^2));
    #newConstraint!(r,obs_con3,:obs_con3);

    # LiDAR connstraint  TODO finish fixing the lidar constraints here NOTE currently not using this
    @NLparameter(mdl, X0_params[j=1:2]==n.X0[j]);

  #  LiDAR_con=@NLconstraint(mdl, [i=1:n.numStatePoints-1], ((r.x[(i+1),1]-X0_params[1])^2+(r.x[(i+1),2]-X0_params[2])^2) <= (c.m.Lr + c.m.L_rd)^2); # not constraining the first state
  #  newConstraint(r,LiDAR_con,:LiDAR_con);
          #   1-3      4

    # constraint on progress on track (no turning around!)
    #=
    if c.t.dir==:posY
      progress_con=@NLconstraint(mdl, [i=1:n.numStatePoints-1], r.x[i,2] <= r.x[(i+1),2]);
      newConstraint!(r,progress_con,:progress_con);
    elseif c.t.dir==:posX
      progress_con=@NLconstraint(mdl, [i=1:n.numStatePoints-1], r.x[i,1] <= r.x[(i+1),1]);
      newConstraint!(r,progress_con,:progress_con);
    end
=#
    # intial optimization
    optimize!(mdl,n,r,s);

    if c.m.activeSafety
      evalConstraints!(n,r);         # setup constraint data
      d=JuMP.NLPEvaluator(mdl);      # initialize feasibility checking functionality
      MathProgBase.initialize(d,[:Grad]);
      checkFeasibility!(mdl,d,n,c,x,pa,r;feas_tol=0.005);
    else
      d=NaN
    end
          #  1      2          3         4
    params=[pa,veh_params, obs_params,X0_params];
    return mdl,n,r,params,x,d
end
"""

checkFeasibility!(mdl,d,n,c,x,params[1],r;feas_tol=0.005);
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 2/13/2017, Last Modified: 4/11/2017 \n
--------------------------------------------------------------------------------------\n
"""
function checkFeasibility!(mdl,d,n,c,x,pa,r;feas_tol::Float64=0.005)

  # simulate vehicle response
  U=zeros(n.numStatePoints,2); # NOTE using numStatePoints so that we can approximate the state at the end point
  SR=copy((x.SA[1]-x.SA[2])/x.dt)*ones(n.numStatePoints,);  # first order approximation of SR
  #JX=0.0;
  U[:,1]=SR;
  #U[:,2]=JX; #NOTE for now assue a constant steering angle
  t0=copy(x.t0); tf=t0+copy(c.m.tp); # tf should be the same as t[end]
  t=t0+r.t_st; # for a fixed final time, this never changes
  sol=simModel(n,pa,x.X0,t,U,t0,tf);

  # pass values to constraint equations
  num=n.numStatePoints*n.numStates + n.numControlPoints*n.numControls; # NOTE we removed t0 as a design variable
  values=zeros(num);
  for st in 1:n.numStates
    for j in 1:n.numStatePoints
     values[linearindex(r.x[j,st])] = sol(t[j])[st];
    end
  end
  #idx=1; #NOTE for now assue a constant steering angle
  #for j in 1:n.numControlPoints;  # assuming JX is zero
  #  values[linearindex(r.u[idx,1])] = SR[idx]; idx+=1;
  #end
  g=zeros(MathProgBase.numconstr(d.m));; # TODO fix -->number of constraints is different from num TODO remove that if we are not doing mpc??
  MathProgBase.eval_g(d,g,values)# is this broken ?
  b=JuMP.constraintbounds(mdl);
  AL=g-b[1]-feas_tol; AU=b[2]+feas_tol-g; # give tolerance

  # check obstacle avoidance constraints
  L1=0;U1=0;  #L2=0;U2=0;  L3=0;U3=0;

  for i in 1:length(r.constraint.name) # note only setup for three obstacles currently
    if r.constraint.name[i]==:obs_con1
      L1=r.constraint.nums[i][end][1];
      U1=r.constraint.nums[i][end][2]-c.m.NF;
    end

#=
    if r.constraint.name[i]==:obs_con2
      L2=r.constraint.nums[i][end][1];
      U2=r.constraint.nums[i][end][2]-c.m.NF;
    end
    if r.constraint.name[i]==:obs_con3
      L3=r.constraint.nums[i][end][1];
      U3=r.constraint.nums[i][end][2]-c.m.NF;
    end
    =#
  end
  #x.AL=[AL[L1:U1];AL[L2:U2];AL[L3:U3]];
#  x.AU=[AU[L1:U1];AU[L2:U2];AU[L3:U3]];
  x.AL=AL[L1:U1];
  x.AU=AU[L1:U1];
  if minimum(sign([x.AL;x.AU]))==-1;x.feasible=false;else;x.feasible=true; end
  nothing
end

"""

--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 4/6/2017, Last Modified: 4/7/2017 \n
--------------------------------------------------------------------------------------\n
"""
function getPlantData!(n,params,x,c,r)
  MsgString = recv(x.s1);L=19;
  MsgIn = zeros(L);
  allidx = find(MsgString->MsgString == 0x20,MsgString);
  allidx = [0;allidx];
  for i in 1:L  #length(allidx)-1
    MsgIn[i] = parse(Float64,String(copy(MsgString[allidx[i]+1:allidx[i+1]-1])));
  end
  x.X0=zeros(n.numStates);
  x.X0[1:7]=MsgIn[1:7];              # global X, global y, lateral speed v, yaw rate r, yaw angle psi
  x.SA=zeros(2);
  x.SA=[copy(x.X0[6]),MsgIn[8]];     # Vehicle steering angle and previous state of steering angle
  x.UX=copy(x.X0[7]);

  N=n.numStates+1;
  x.Xobs=MsgIn[N:N+x.numObs-1];               # ObsX
  x.Yobs=MsgIn[N+x.numObs:N+2*x.numObs-1];    # ObsY
  x.A= MsgIn[N+2*x.numObs:N+3*x.numObs-1];    # ObsR
  x.B=x.A;
  x.runJulia=MsgIn[N+3*x.numObs];
  x.t0=MsgIn[N+3*x.numObs+1];
  n.mpc.t0_actual=copy(x.t0);

  # update drivers steering rate
  SR=copy((x.SA[1]-x.SA[2])/x.dt);
  if SR>0.086;SR=0.086;elseif SR<-0.086;SR=0.086;end # saturate SR
  setvalue(params[2][3],SR)

  # update obstacle feild
#  for i in 1:length(x.A)
  i=1;
  setvalue(params[3][1][i],x.A[i]);
  setvalue(params[3][2][i],x.B[i]);
  setvalue(params[3][3][i],x.Xobs[i]);
  setvalue(params[3][4][i],x.Yobs[i]);
#  end

  # update case for plotting
  if r.eval_num==1; c.o.A=[];c.o.B=[];c.o.X0=[];c.o.Y0=[]end # reset

  if isempty(c.o.X0) # update obstacle data
    push!(c.o.A,copy(x.A[1]));
    push!(c.o.B,copy(x.B[1]));
    push!(c.o.X0,copy(x.Xobs[1]));
    push!(c.o.Y0,copy(x.Yobs[1]));
    push!(c.o.s_y,0.0);
    push!(c.o.s_x,0.0);
  elseif  (x.Xobs!=c.o.X0[end]) && (x.Yobs!=c.o.Y0[end])
    push!(c.o.A,copy(x.A[1]));
    push!(c.o.B,copy(x.B[1]));
    push!(c.o.X0,copy(x.Xobs[1]));
    push!(c.o.Y0,copy(x.Yobs[1]));
    push!(c.o.s_y,0.0);
    push!(c.o.s_x,0.0);
  end

  x.SA_first=zeros(n.numStatePoints,);
  for jj in 1:n.numStatePoints
    x.SA_first[jj]=copy(x.SA[1]) + jj*copy(x.SA[1]-x.SA[2]);  # first order approximation of SA
  end

  if c.m.followDriver # first order
    for jj in 1:n.numStatePoints
      setvalue(params[2][2][jj], x.SA_first[jj])
    end
  end
  nothing
end

"""

--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 4/6/2017, Last Modified: 4/7/2017 \n
--------------------------------------------------------------------------------------\n
"""
function sendOptData!(n,r,x)
  if !x.feasible || (x.infeasible_counter < x.infeasible_counter_max)
    println("Sending Optimized Steering Angle Commands to MATLAB");
    if r.dfs_opt[r.eval_num-1][:status][1]==:Optimal # if infeasible -> let user control -> the 3 is a flag in Matlab that tells the system not to use signals from Autonomy
      MsgOut = [x.sa_sample;r.dfs_opt[r.eval_num-1][:t_solve];2;1]
    else
      MsgOut = [x.sa_sample;r.dfs_opt[r.eval_num-1][:t_solve];3;1];
    end
    # send UDP packets to client side
    TJ=copy(n.mpc.t0_actual)
    MsgOut = [MsgOut;Float64(r.eval_num-1);TJ];
  else
    MsgOut=[copy(x.SA[1])*ones(31);0.01;1;1]
  end
  MsgOutString = ' ';
  for j in 1:length(MsgOut)
      MsgOutString = string(MsgOutString,' ',MsgOut[j]);
  end
  MsgOutString = string(MsgOutString," \n");
  send(x.s2,ip"141.212.141.245",36881,MsgOutString);  # change this to the ip where you are running Simulink!
  nothing
end
"""

--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 1/27/2017, Last Modified: 4/7/2017 \n
--------------------------------------------------------------------------------------\n
"""
function sharedControl!(mdl,n,r,s,params,x)
  # rerun optimization
  status=autonomousControl!(mdl,n,r,s,params[1]);

  # sample solution
  sp_SA=Linear_Spline(r.t_st,r.X[:,6]);
  t_sample=Ranges.linspace(0.0,3.0,31);
  x.sa_sample=sp_SA[t_sample];

  # update status for Matlab
  if status!=:Infeasible; x.status=1.; else x.status=0.; end
  nothing
end

"""
evalNum()
# to extract data from a particular iteration number
# will not work unless the data was saved; i.e., s.save = true
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 2/21/2017, Last Modified: 2/21/2017 \n
--------------------------------------------------------------------------------------\n
"""
function evalNum(Idx)
  eval_num=0;
  for i in 1:length(r.dfs_opt)
    if r.dfs_opt[i][:iter_num][1]==Idx
      eval_num=i;
      break
    end
  end
  print(eval_num);
  s=Settings(;format=:png,MPC=false);
  cd("results/test1"); allPlots(n,r,s,eval_num); cd(main_dir)
end

end # module


text:
using CompScienceMeshes
using BEAST

m = meshsphere(1.0, 0.45)
# m = CompScienceMeshes.rotate(m, rand(3))
X = raviartthomas(m)

Y = buffachristiansen(m)
# Y = X
Y = BEAST.subset(Y,1:135)
geo = geometry(Y)

els, ad = BEAST.assemblydata(Y)
Base.length(it::BEAST.ADIterator) = something(findfirst(x->x[1]<1, it.ad.data[:,it.r,it.c]), size(it.ad.data,1)+1)-1

Ad = [Int[] for c in 1:length(els)]
for c in 1:length(els)
    for r in 1:3
        length(ad[c,r]) == 0 && continue
        push!(Ad[c], maximum(m for (m,a) in ad[c,r]))
    end
end
top = maximum.(Ad)
bot = minimum.(Ad)

using Plots

plot()
plot!(top)
plot!(bot)

import PlotlyJS

t1 = patch(Mesh(geo.mesh.vertices, geo.mesh.faces[1:150]), rand(150))
P = [p[i] for p in Y.pos, i in 1:3]
t2 = PlotlyJS.scatter3d(x=P[1:38,1], y=P[1:38,2], z=P[1:38,3])
t3 = PlotlyJS.scatter3d(x=P[264:264,1], y=P[264:264,2], z=P[264:264,3])
t4 = PlotlyJS.scatter3d(x=P[1:180,1], y=P[1:180,2], z=P[1:180,3])
PlotlyJS.plot(t4)

PlotlyJS.plot([t1,t2,t3])

ext: jl


text:
using Test

@testset "Solution.jl" begin
    using Sum_of_Powers: Solution, err
    s = Solution(15,16,(14,9,8,6,3))
    @test string(s) == "15^16=>{3,6,8,9,14}"
    @test err(s) == BigInt(4388317701585002815)
    @test err(Solution(24,16,Tuple(1:23))) == BigInt(400030525071869538932)
end

@testset "OnePositions.jl" begin
    using Sum_of_Powers: OnePositions
    test_one_positions(x::T) where {T<:Integer} = sum(map(x_->T(2)^(x_-1),OnePositions(x))) == x
    @test test_one_positions(0)
    @test test_one_positions(1)
    @test test_one_positions(2)
    @test test_one_positions(Int16(2)^15-1)
    @test test_one_positions(Int32(2)^31-1)
    @test test_one_positions(Int64(2)^63-1)
    @test test_one_positions(Int128(2)^127-1)
    for x in rand(1:Int64(2)^63-1,10)
        @test test_one_positions(x)
    end
    for x in rand(1:Int128(2)^127-1,10)
        @test test_one_positions(x)
    end
    a = [1,2,3,4,5,6,7]
    @test setdiff(a[OnePositions(10)],[2,4])==[] # getindex
    b = falses(8)
    b[OnePositions(10)] = 1
    @test b==[0,1,0,1,0,0,0,0] # setindex!
end

@testset "track_best.jl" begin
    using Sum_of_Powers: Tracker, is_error_zero
    s = 10
    n = 5
    lhs = BigInt(s)^n
    a = 0
    best = Tracker{Int}(a,lhs)
    best(31,-lhs-1)
    (best_a,best_e) = best()
    @test best_a==a
    @test ~is_error_zero(best)
    new_a = 255
    best(new_a,BigInt(0))
    (best_a,best_e) = best()
    @test best_a==new_a
    @test is_error_zero(best)
end

@testset "search.jl" begin
    using Sum_of_Powers: search, Best, SubSet, SlidingWindow, MabeyBest
    a, e = search(Best(),20,10)
    @test e == BigInt(8584238000)
    @test a == [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 18, 19]
    search(Best(),20,10,[1,2,3,4]) # should ignore array

    ss = SubSet([x for x in 1:8])
    s = 20
    n = 16
    a = [x for x in 1:19]
    known_e = 168716336388812100534 # for {19,18,17,16,15,14,13,12,11,10,9,8,7,6,5,4,3,2,1}
    a,e = search(ss,s,n,a)
    @test e==known_e
    @test setdiff(a,[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19])==[]

    ss = SubSet([14,13,4,5,6,10])
    s = 16
    n = 10
    a = [15,12,11]
    known_e = -2106698499
    a,e = search(ss,s,n,a)
    @test e==known_e
    @test setdiff(a,[15,14,13,12,11,10])==[]
    sol = Solution(s,n,a,false)
    @test err(sol) == known_e
    # best 28^16=>{27,26,24,23,20,19,18,17},e=-3923372792424650116
    s = 28
    n = 16
    a = Vector{Int}()
    sw = SlidingWindow(10)
    a,e = search(sw,s,n,a)
    @test e==-3923372792424650116 # SlidingWindow found the best!
    @test setdiff(a,[27, 26, 24, 23, 20, 19, 18, 17])==[]

    mb = MabeyBest()
    # 2^3=>{1},e=7
    s = 2
    n = 3
    a,e = search(mb,s,n)
    @test e==7
    @test issetequal(a,[1])
    # 5^3=>{4,3,2,1},e=25
    s = 5
    n = 3
    a,e = search(mb,s,n)
    @test e==25
    @test issetequal(a,[4,3,2,1])
    # 6^3=>{5,4,3},e=0
    s = 6
    n = 3
    a,e = search(mb,s,n)
    @test e==0
    @test issetequal(a,[5,4,3])
    # 11^6=>{10,9,8},e=-22024
    s = 11
    n = 6
    a,e = search(mb,s,n)
    @test e==-22024
    @test issetequal(a,[10,9,8])
    # 25^6=>{23,18,17,16,15,13,12,10,9,8,7,6,5,3,2,1},e=0
    s = 25
    n = 6
    a,e = search(mb,s,n)
    @test e==0
    @test issetequal(a,[23,18,17,16,15,13,12,10,9,8,7,6,5,3,2,1])
    # 28^40=>{27,26,25,24,23,22,21,20,19,18,17,16,15,14,13,12,11,10,9,8,7,6,5,4,3,2,1},e=5400443268136291573096133732738005254949693470241623252754
    s = 28
    n = 40
    a,e = search(mb,s,n)
    @test e==BigInt(5400443268136291573096133732738005254949693470241623252754)
    @test issetequal(a,[27,26,25,24,23,22,21,20,19,18,17,16,15,14,13,12,11,10,9,8,7,6,5,4,3,2,1])
    # 28^16=>{27,26,24,23,20,19,18,17},e=-3923372792424650116
    s = 28
    n = 16
    a,e = search(mb,s,n)
    @test e==BigInt(-3923372792424650116)
    @test issetequal(a,[27,26,24,23,20,19,18,17])
    #25^16=>{24,23,22,21,20},e=-63411027820192788929
    s = 25
    n = 16
    a,e = search(mb,s,n)
    @test e==BigInt(-63411027820192788929)
    @test issetequal(a,[24,23,22,21,20])
end


text:
module MaxHelpingHandGradientDustTrigger

using ..Ahorn, Maple

@mapdef Trigger "MaxHelpingHand/GradientDustTrigger" GradientDustTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight,
    gradientImage::String="MaxHelpingHand/gradientdustbunnies/bluegreen", scrollSpeed::Number=50.0)

const placements = Ahorn.PlacementDict(
    "Gradient Dust (max480's Helping Hand)" => Ahorn.EntityPlacement(
        GradientDustTrigger,
        "rectangle"
    )
)

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
module QuantumESPRESSOExpress

using AbInitioSoftwareBase: AbInitioSoftware

import Express: current_software

export QE

struct QuantumESPRESSO <: AbInitioSoftware end
const QE = QuantumESPRESSO

currentsoftware() = QE()

include("ConvergenceTestWorkflow/ConvergenceTestWorkflow.jl")
include("EquationOfStateWorkflow/EquationOfStateWorkflow.jl")
include("PhononWorkflow/PhononWorkflow.jl")

end


text:
using Compose, Viznet
using Test

@testset "sq lattice" begin
    lt = rectlattice(10, 20, 1.0, 1.5)
    @test size(lt, 1) == 10
    @test bravais_size(lt) == (10, 20)
    @test size(lt, 2) == 20
    @test size(lt) == (10, 20)
    @test unit(lt) ≈ 1/30
    @test vec_a(lt) ≈ [1.0, 0.0]
    @test vec_b(lt) ≈ [0.0, 1.5]
end

@testset "bonds" begin
    lt = rectlattice(4, 5, 1.0, 2.0)
    b1s = bonds(lt)
    @test length(b1s) == 31
    ti = true
    for b in b1s
        ti = ti && isconnected(lt, b...)
    end
    @test ti
end


text:
using vOptSpecific

id = load2OSP( "2osp04.dat" )
solver = OSP_VanWassenhove1980( )
z1, z2 , S = vSolve( id , solver )

max_stars_repo_name: vOptSolver/vOptSpecific

max_issues_repo_name: vOptSolver/vOptSpecific

max_forks_repo_name: vOptSolver/vOptSpecific


text:
module Step
    export oper_array
    export gcdi
    export lcmu
    export som
    export maxi
    export mini

    gcdi = gcd
    lcmu = lcm
    som = +
    maxi = max
    mini = min

    oper_array(fct, arr, init) = map( x -> init=fct(init,x), arr )

end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
@info "Testing CustomDataset..."
@testset "CustomDataset" begin
    x = rand(100, 8)
    getdata(idx) = x[idx, :]
    dset = CustomDataset(size(x, 1), getdata)

    @test size(dset) == (size(x, 1),)
    @test length(dset) == size(x, 1)
    @test firstindex(dset) == 1
    @test lastindex(dset) == size(x, 1)
    @test lastindex(dset) == size(x, 1)
    @test dset[1] == x[1, :]
    @test dset[3:8] == [x[i, :] for i = 3:8]
    @test dset[:] == [x[i, :] for i = 1:size(x, 1)]
end


text:
using SolidStateDetectors
using IntervalSets
using SpecialFunctions
T=Float32

struct CustomChargeDensity{T, N} <: SolidStateDetectors.AbstractChargeDensity{T}
	G1::SolidStateDetectors.AbstractGeometry{T,N}
	G2::SolidStateDetectors.AbstractGeometry{T,N}
	L1::SolidStateDetectors.AbstractChargeDensity{T}
	L2::SolidStateDetectors.AbstractChargeDensity{T}
	G3::SolidStateDetectors.AbstractGeometry{T,N}
	L3::SolidStateDetectors.AbstractChargeDensity{T}
	G4::SolidStateDetectors.AbstractGeometry{T,N}
	L4::SolidStateDetectors.AbstractChargeDensity{T}
	G5::SolidStateDetectors.AbstractGeometry{T,N}
	L5::SolidStateDetectors.AbstractChargeDensity{T}
	G6::SolidStateDetectors.AbstractGeometry{T,N}
	L6::SolidStateDetectors.AbstractChargeDensity{T}
	G7::SolidStateDetectors.AbstractGeometry{T,N}
	L7::SolidStateDetectors.AbstractChargeDensity{T}
	G8::SolidStateDetectors.AbstractGeometry{T,N}
	L8::SolidStateDetectors.AbstractChargeDensity{T}
	G9::SolidStateDetectors.AbstractGeometry{T,N}
	L9::SolidStateDetectors.AbstractChargeDensity{T}
end

function SolidStateDetectors.get_charge_density(cdm::CustomChargeDensity{T}, pt::SolidStateDetectors.AbstractCoordinatePoint{T})::T where {T}
	(pt in cdm.G1) * SolidStateDetectors.get_charge_density(cdm.L1, pt) +
	(pt in cdm.G2) * SolidStateDetectors.get_charge_density(cdm.L2, pt) +
	(pt in cdm.G3) * SolidStateDetectors.get_charge_density(cdm.L3, pt) +
	(pt in cdm.G4) * SolidStateDetectors.get_charge_density(cdm.L4, pt) +
	(pt in cdm.G5) * SolidStateDetectors.get_charge_density(cdm.L5, pt) +
	(pt in cdm.G6) * SolidStateDetectors.get_charge_density(cdm.L6, pt) +
	(pt in cdm.G7) * SolidStateDetectors.get_charge_density(cdm.L7, pt) +
	(pt in cdm.G8) * SolidStateDetectors.get_charge_density(cdm.L8, pt) +
	(pt in cdm.G9) * SolidStateDetectors.get_charge_density(cdm.L9, pt)
end

G1 = SolidStateDetectors.Tube{T}(0..0.01627006165897972,0..2π,0..0.002,CartesianPoint{T}(-1.2390881971262907e-19, 9.912705577010326e-19, 0))
G2 = SolidStateDetectors.Tube{T}(0..0.01627006165897972,0..2π,0..1e-05,CartesianPoint{T}(-1.2390881971262907e-19, 9.912705577010326e-19, 0))
CD1 = SolidStateDetectors.CylindricalChargeDensity{T}((0.0f0,0.0f0,50000000000.0),(60000000000.0,0.0f0,0.0f0))
CD2 = SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,-1.9e+16),(0.0f0,0.0f0,0.0f0))
G3 = SolidStateDetectors.HexagonalPrism{T}(0, 0.00515, 2e-05,CartesianPoint{T}(0, 0, 0.00199),0)
CD3= SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,1.9e+16),(0.0f0,0.0f0,0.0f0))
G4 = SolidStateDetectors.HexagonalPrism{T}(0, 0.00515, 2e-05,CartesianPoint{T}(0.009020061658979717, 0.0, 0.00199),0)
CD4= SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,1.9e+16),(0.0f0,0.0f0,0.0f0))
G5 = SolidStateDetectors.HexagonalPrism{T}(0, 0.00515, 2e-05,CartesianPoint{T}(0.0045100308294898596, 0.0078116025403784425, 0.00199),0)
CD5= SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,1.9e+16),(0.0f0,0.0f0,0.0f0))
G6 = SolidStateDetectors.HexagonalPrism{T}(0, 0.00515, 2e-05,CartesianPoint{T}(-0.004510030829489857, 0.007811602540378445, 0.00199),0)
CD6= SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,1.9e+16),(0.0f0,0.0f0,0.0f0))
G7 = SolidStateDetectors.HexagonalPrism{T}(0, 0.00515, 2e-05,CartesianPoint{T}(-0.009020061658979717, 1.1046389638781276e-18, 0.00199),0)
CD7= SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,1.9e+16),(0.0f0,0.0f0,0.0f0))
G8 = SolidStateDetectors.HexagonalPrism{T}(0, 0.00515, 2e-05,CartesianPoint{T}(-0.004510030829489863, -0.007811602540378441, 0.00199),0)
CD8= SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,1.9e+16),(0.0f0,0.0f0,0.0f0))
G9 = SolidStateDetectors.HexagonalPrism{T}(0, 0.00515, 2e-05,CartesianPoint{T}(0.004510030829489859, -0.0078116025403784425, 0.00199),0)
CD9= SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,1.9e+16),(0.0f0,0.0f0,0.0f0))

ccdm = CustomChargeDensity(G1,G2,CD1,CD2,G3,CD3,G4,CD4,G5,CD5,G6,CD6,G7,CD7,G8,CD8,G9,CD9)

max_stars_repo_name: leenderthayen/NDP

max_issues_repo_name: leenderthayen/NDP

max_forks_repo_name: leenderthayen/NDP


text:
using Batman, Test

m = CountingExperiment()
add_component!(m, "Signal", 20.0)
add_component!(m, "Bkg 1", 30.0; σ = 0.1)
add_component!(m, "Bkg 2", 30.0; σ = 45.0)
add_component!(m, "Bkg 3", 12.0; σ = 12.0)

set_counts!(m, 100)

results = minimize!(m)

compute_profiled_uncertainties!(results; σ=1)

signal = getparam(m, "Signal")

@test signal.fit ≈ 18 rtol=0.5


text:
"""
Before performing disk IO, data is manipulated within a `Pgae`.
Data is broken up into `Block` objects, which are then saved to file.
"""
module Pages
  using HelloDB.FileMgrs
  using HelloDB.Blocks
  using HelloDB.Data: INT_SIZE, CHAR_SIZE, BLOCK_SIZE
  using HelloDB: DB

  ##############################################################################
  # Exports
  ##############################################################################

  export Page
  export string_size, getint, setint, getstring, setstring, append

  ##############################################################################
  # Implementation
  ##############################################################################

  type Page
    contents::IOBuffer
    filemgr::FileMgr

    """
    Page() is initialized with the Database file manage (`DB.filemgr`).
    For testing purposes, one may initialize a `Page` with an empty
    filemanager which will not write to disk by negating the
    `realfile` initialization parameter: `Page(false)`.
    """
    function Page(realfile=true)
      buffer = IOBuffer(BLOCK_SIZE)

      if realfile
        filemgr = DB.filemgr
      else
        filemgr = FileMgr()
      end

      new(buffer, filemgr)
    end
  end

  """
  Append `page` contents to file `filename`.
  """
  function append(page::Page, filename)
    FileMgrs.append(page.filemgr, filename, page.contents)
  end

  """
  Gets the first byte in `page.contents` buffer as `Int`.
  (This may represnt the length of string starting at the pointer position.)
  """
  function getint(page::Page, offset::Int)
    seek(page.contents, offset)
    read(page.contents, Int)
  end

  """
  Like `getint()` but for a string of `Char`.
  First we get the string's expected length (represented by an int at the
  beginning of the char sequence), next we red the chars into a string.
  """
  function getstring(page::Page, offset::Int)
    seek(page.contents, offset)
    str_len = read(page.contents, Int)
    # We assume all chars are stored on disk as single-byte, then convert these
    # hex values to an array of chars in memory, then to a string
    # (grep TODO#charsize).
    AbstractString(map(x->Char(x), readbytes(page.contents, str_len)))
  end

  """
  Fill contents of `page` based on `block` filename and position.
  """
  function Base.read(page::Page, block::Block)
    read(page.filemgr, block, page.contents)
  end

  """
  Insert `val` into `page.contents` byte buffer at offset `offset`.
  """
  function setint(page::Page, offset::Int, val::Int)
    seek(page.contents, offset)
    write(page.contents, val)
  end

  """
  Like `setint()` but for a string of `Char`.
  """
  function setstring(page::Page, offset::Int, str::AbstractString)
    seek(page.contents, offset)
    write(page.contents, length(str))
    write(page.contents, IOBuffer(str))
  end

  """
  Strings will be represented as arrays of the form `[str_length, chars...]`.
  """
  function string_size(n)
    INT_SIZE + (n * CHAR_SIZE)
  end

  """
  Write to `block` from `page` contents.
  """
  function Base.write(page::Page, block::Block)
    write(page.filemgr, block, page.contents)
  end
end


text:
module MatrixModule
function readMatrix(filename::String,dimensions::Tuple{Integer,Integer}=(0,0))
  floats = reinterpret(Float64,read(filename))

  if dimensions==(0,0) # try quadratic matrix
    dim = Int64(sqrt(length(floats)))
    if dim*dim != length(floats)
      error("Non-square matrix. Specify dimensions.")
    end
    dimensions = (dim,dim)
  end

  return reshape(floats,dimensions)
end

function ASCIIArt(matrix::Matrix)
  result = Matrix{Char}(size(matrix))
  for i in 1:length(matrix)
    if matrix[i] >= 10e-0
      result[i] = 'M'
    elseif matrix[i] >= 10e-1
      result[i] = 'N'
    elseif matrix[i] >= 10e-2
      result[i] = '@'
    elseif matrix[i] >= 10e-3
      result[i] = '%'
    elseif matrix[i] >= 10e-4
      result[i] = '#'
    elseif matrix[i] >= 10e-5
      result[i] = 't'
    elseif matrix[i] >= 10e-6
      result[i] = '+'
    elseif matrix[i] >= 10e-7
      result[i] = ';'
    elseif matrix[i] >= 10e-8
      result[i] = ':'
    elseif matrix[i] >= 10e-9
      result[i] = ','
    elseif matrix[i] >= 10e-10
      result[i] = '.'
    else
      result[i] = ' '
    end
  end
  return result
end

function writePixmap(filename::String,matrix::Matrix)
  # taking the specification from wikipedia
  dimensions = size(matrix)
  colors = ['M' => "#ff0000",
            'N' => "#ff3300",
            '@' => "#ff3333",
            '%' => "#ff6633",
            '#' => "#ff6666",
            't' => "#ff9966",
            '+' => "#ff9999",
            ';' => "#ffbb99",
            ':' => "#ffbbbb",
            ',' => "#ffeebb",
            '.' => "#ffeeee",
            ' ' => "#ffffff",
	    '0' => "#bbbbbb", # zero by matrix format (blocking)
	   ]
  asciiart = ASCIIArt(matrix)
  open(filename,"w") do fn
    write(fn,"/* XPM */\n")
    write(fn,"static char* $(split(filename,'.')[1])[] = {\n")
    write(fn,"\"$(dimensions[2]) $(dimensions[1]) $(length(colors)) 1\",\n")
    for (sym,color) in colors
      write(fn,"\"$sym c $color\",\n")
    end
    for line in 1:dimensions[1]
      write(fn,"\"$(join(asciiart[line,:]))\",\n")
    end
  end
end
end


text:
import NumDiffTools
norm(x) = maximum(abs.(x))
using Test
using ForwardDiff

@info "Checking derivative(f,x) of scalar function"
err = NumDiffTools.derivative(sin, 2.0) - cos(2.0)
@test norm(err) < 1e-12

@info "Checking derivative(f,x) of vector function"
f(x) = [sin(x), cos(x)]
f!(y,x) = begin
    y[1] = sin(x)
    y[2] = cos(x)
    nothing
end
df(x) = [cos(x), -sin(x)]
err = NumDiffTools.derivative(f, 2.0) - df(2.0)
@test norm(err) < 1e-12

@info "Checking derivative(f!, y, x)"
y = zeros(2)
err =  NumDiffTools.derivative(f!, y, 2.0) - df(2.0)
@test norm(err) < 1e-12

@info "Checking derivative!(result, f, x)"
result = zeros(2)
NumDiffTools.derivative!(result, f, 2.0)
err = result - df(2)
@test norm(err) < 1e-12

@info "Checking derivative!(result, f, y, x)"
NumDiffTools.derivative!(result,f!, y, 2.0)
err = result - df(2.0)
@test norm(err) < 1e-12

@info "Checking gradient(f,x)"
f(x) = prod(cos.(x).^2 - sin.(x).^2)
gradf(x) = [ - 4*cos(x[1])*sin(x[1])*(cos(x[2])^2 - sin(x[2])^2),
             - 4*cos(x[2])*sin(x[2])*(cos(x[1])^2 - sin(x[1])^2)
             ]
err = NumDiffTools.gradient(f, [1.5, 2.5]) - gradf([1.5, 2.5])
@test norm(err) < 1e-12

@info "Checking gradient!(result, f, x)"
NumDiffTools.gradient!(result, f, [1.5, 2.5])
err = result - gradf([1.5, 2.5])
@test norm(err) < 1e-12

@info "checking jacobian(f, x)"
f(x) = [ prod(cos.(x).^2 - sin.(x).^2),
         sum(cos.(x).^2 - sin.(x).^2)]
f!(y, x) = begin
    y .= f(x)
    nothing
end
J(x) = [-4*cos(x[1])*sin(x[1])*(cos(x[2])^2 - sin(x[2])^2)  -4*cos(x[2])*sin(x[2])*(cos(x[1])^2 - sin(x[1])^2);
        (-2*cos(x[1])*sin(x[1]) - 2*sin(x[1])*cos(x[1]))  (-2*cos(x[2])*sin(x[2]) - 2*sin(x[2])*cos(x[2]))
        ]
err = NumDiffTools.jacobian(f, [1.5, 2.5]) - J([1.5, 2.5])
@test norm(err) < 1e-12

@info "checking jacobian(f!, y, x)"
err = NumDiffTools.jacobian(f!, y, [1.5, 2.5]) - J([1.5, 2.5])
@test norm(err) < 1e-12

@info "checking jacobian!(result, f, x)"
result = zeros(2,2)
NumDiffTools.jacobian!(result, f, [1.5, 2.5])
err = result - J([1.5, 2.5])
@test norm(err) < 1e-12

@info "Checking jacobian!(result, f!, y, x)"
result = zeros(2,2)
NumDiffTools.jacobian!(result, f!, y, [1.5, 2.5])
err = result - J([1.5, 2.5])
@test norm(err) < 1e-12

@info "Checking hessian(f, x) for scalar functions of real numbers"
err = NumDiffTools.hessian(sin, 1.0) + sin(1.0)
@test norm(err) < 1e-10

@info "checkking hessian(f, x) for scalar functions of vectors"
f(x) = prod(cos.(x).^2 - sin.(x).^2)

err = NumDiffTools.hessian(f, [1.5, 2.5]) - ForwardDiff.hessian(f, [1.5, 2.5])
@test norm(err) < 1e-10

@info "Checking hessian!(result, f, x)"
result = zeros(2,2)
NumDiffTools.hessian!(result, f, [1.5, 2.5])
err = result - ForwardDiff.hessian(f, [1.5, 2.5])
@test norm(err) < 1e-10


text:
include("test_read_write.jl")
include("test_split_trigger.jl")
include("test_split_time1.jl")
include("test_split_time2.jl")
include("test_eq.jl")
include("test_channel_select.jl")
include("test_decode_status_channel.jl")


text:
Base.strides(A::AxisArray) = strides(parent(A))

"""
    AxisVector

A vector whose indices have keys.

## Examples
```jldoctest
julia> using AxisIndices

julia> AxisVector([1, 2], [:a, :b])
2-element AxisArray(::Vector{Int64}
  • axes:
     1 = [:a, :b]
)
      1
  :a  1
  :b  2

```
"""
const AxisVector{T,P<:AbstractVector{T},Ax} = AxisArray{T,1,P,Tuple{Ax}}

function AxisVector{T}(x::AbstractVector{T}, ks::AbstractVector) where {T}
    axis = Axis(ks, axes(x, 1))
    return AxisArray{T,1,typeof(x),Tuple{typeof(axis)}}(x, (axis,))
end

AxisVector(x::AbstractVector{T}, ks::AbstractVector) where {T} = AxisVector{T}(x, ks)

AxisVector(x::AbstractVector) = AxisArray(x)

function AxisVector{T}() where {T}
    return AxisArray{T,1,Vector{T},Tuple{SimpleAxis{Int,OneToMRange{Int}}}}(
        T[], (SimpleAxis(OneToMRange(0)),)
    )
end

@inline function append_axis!(x::Axis{K,I,Ks,Inds}, y) where {K,I,Ks,Inds}
    if Ks <: AbstractRange
        set_length!(values(x), length(x) + length(y))
    else
        if any(in(keys(x)), keys(y))
            error("Cannot append axis keys that are not unique from each other.")
        else
            append!(keys(x), keys(y))
        end
    end
end
function Base.append!(A::AxisVector{T,V,Ax}, collection) where {T,V,Ax}
    if Ax <: Axis
        append_axis!(axes(A, 1), axes(collection, 1))
    else
        set_length!(axes(A, 1), length(A) + length(collection))
    end
    append!(parent(A), collection)
    return A
end

function Base.pop!(A::AxisVector)
    shrink_last!(axes(A, 1), 1)
    return pop!(parent(A))
end

function Base.popfirst!(A::AxisVector)
    popfirst_axis!(axes(A, 1))
    return popfirst!(parent(A))
end

function Base.reverse(x::AxisVector)
    p = reverse(parent(x))
    return initialize_axis_array(p, (reverse_keys(axes(x, 1), axes(p, 1)),))
end

"""
    deleteat!(a::AxisVector, arg)

Remove the items corresponding to `A[arg]`, and return the modified `a`. Subsequent
items are shifted to fill the resulting gap. If the axis of `a` is an `SimpleAxis`
then it is shortened to match the length of `a`.

## Examples
```jldoctest
julia> using AxisIndices

julia> x = AxisArray([1, 2, 3, 4]);

julia> deleteat!(x, 3)
3-element AxisArray(::Vector{Int64}
  • axes:
     1 = 1:3
)
     1
  1  1
  2  2
  3  4

julia> x = AxisArray([1, 2, 3, 4], ["a", "b", "c", "d"]);

julia> keys.(axes(deleteat!(x, "c")))
(["a", "b", "d"],)

```
"""
function Base.deleteat!(A::AxisVector{T,P,Ax}, arg) where {T,P,Ax}
    if Ax<:Axis
        inds = to_index(axes(A, 1), arg)
        deleteat!(keys(axes(A, 1)), inds)
        shrink_last!(parent(axes(A, 1)), length(inds))
        deleteat!(parent(A), inds)
        return A
    else
        inds = to_index(axes(A, 1), arg)
        shrink_last!(axes(A, 1), length(inds))
        deleteat!(parent(A), inds)
        return A
    end
end

function Base.insert!(A::AxisVector, index, item)
    if can_change_size(A)
        axis = axes(A, 1)
        unsafe_insert!(parent(A), axis, to_index(axis, index), item)
        return A
    else
        throw(MethodError(insert!, (A, index, item)))
    end
end

function unsafe_insert!(data::AbstractVector{T}, axis, index::Int, item::I) where {T,I}
    unsafe_insert!(data, axis, index, convert(T, item))
    return nothing
end

function unsafe_insert!(data::AbstractVector{T}, axis, index::Int, item::I) where {T,I<:T}
    grow_last!(axis, 1)
    insert!(data, index, item)
    return nothing
end

function Base.resize!(x::AxisVector, n::Integer)
    resize!(parent(x), n)
    resize_last!(axes(x, 1), n)
    return x
end

###
### Matrix methods
###
const AxisMatrix{T,P<:AbstractMatrix{T},Ax1,Ax2} = AxisArray{T,2,P,Tuple{Ax1,Ax2}}

"""
    rot180(A::AxisMatrix)

Rotate `A` 180 degrees, along with its axes keys.

## Examples
```jldoctest
julia> using AxisIndices

julia> a = AxisArray([1 2; 3 4], ["a", "b"], ["one", "two"]);

julia> b = rot180(a);

julia> keys.(axes(b))
(["b", "a"], ["two", "one"])

julia> c = rotr90(rotr90(a));

julia> keys.(axes(c))
(["b", "a"], ["two", "one"])

julia> a["a", "one"] == b["a", "one"] == c["a", "one"]
true
```
"""
function Base.rot180(x::AxisMatrix)
    p = rot180(parent(x))
    axs = (reverse_keys(axes(x, 1), axes(p, 1)), reverse_keys(axes(x, 2), axes(p, 2)))
    return AxisArray{eltype(p),2,typeof(p),typeof(axs)}(p, axs)
end

"""
    rotr90(A::AxisMatrix)

Rotate `A` right 90 degrees, along with its axes keys.

## Examples
```jldoctest
julia> using AxisIndices

julia> a = AxisArray([1 2; 3 4], ["a", "b"], ["one", "two"]);

julia> b = rotr90(a);

julia> keys.(axes(b))
(["one", "two"], ["b", "a"])

julia> a["a", "one"] == b["one", "a"]
true
```
"""
function Base.rotr90(x::AxisMatrix)
    p = rotr90(parent(x))
    axs = (assign_indices(axes(x, 2), axes(p, 1)), reverse_keys(axes(x, 1), axes(p, 2)))
    return AxisArray{eltype(p),2,typeof(p),typeof(axs)}(p, axs)
end

"""
    rotl90(A::AxisMatrix)

Rotate `A` left 90 degrees, along with its axes keys.

## Examples
```jldoctest
julia> using AxisIndices

julia> a = AxisArray([1 2; 3 4], ["a", "b"], ["one", "two"]);

julia> b = rotl90(a);

julia> keys.(axes(b))
(["two", "one"], ["a", "b"])

julia> a["a", "one"] == b["one", "a"]
true

```
"""
function Base.rotl90(x::AxisMatrix)
    p = rotl90(parent(x))
    axs = (reverse_keys(axes(x, 2), axes(p, 1)), assign_indices(axes(x, 1), axes(p, 2)))
    return AxisArray{eltype(p),2,typeof(p),typeof(axs)}(p, axs)
end

###
### VecOrMat
###
const AxisVecOrMat{T} = Union{<:AxisMatrix{T},<:AxisVector{T}}

###
### reduce
###
function reconstruct_reduction(old_array, new_array, dims)
    return AxisArray(new_array, reduce_axes(axes(old_array), axes(new_array), dims))
end
reconstruct_reduction(old_array, new_array, dims::Colon) = new_array

function Base.mapreduce(f1, f2, a::AxisArray; dims=:, kwargs...)
    return reconstruct_reduction(a, Base.mapreduce(f1, f2, parent(a); dims=dims, kwargs...), dims)
end

function Base.extrema(A::AxisArray; dims=:, kwargs...)
    return reconstruct_reduction(A, Base.extrema(parent(A); dims=dims, kwargs...), dims)
end

if VERSION > v"1.2"
    function Base.has_fast_linear_indexing(x::AxisArray)
        return Base.has_fast_linear_indexing(parent(x))
    end
end

for f in (:mean, :std, :var, :median)
    @eval function Statistics.$f(a::AxisArray; dims=:, kwargs...)
        return reconstruct_reduction(a, Statistics.$f(parent(a); dims=dims, kwargs...), dims)
    end
end

"""
    reshape(A::AxisArray, shape)

Reshape the array and axes of `A`.

## Examples
```jldoctest
julia> using AxisIndices

julia> A = reshape(AxisArray(Vector(1:8), [:a, :b, :c, :d, :e, :f, :g, :h]), 4, 2);

julia> axes(A)
(Axis([:a, :b, :c, :d] => SimpleAxis(1:4)), SimpleAxis(1:2))

julia> axes(reshape(A, 2, :))
(Axis([:a, :b] => SimpleAxis(1:2)), SimpleAxis(1:4))

```
"""
function Base.reshape(A::AxisArray, shp::NTuple{N,Int}) where {N}
    p = reshape(parent(A), shp)
    return AxisArray(p, reshape_axes(naxes(A, Val(N)), axes(p)))
end

function Base.reshape(A::AxisArray, shp::Tuple{Vararg{Union{Int,Colon},N}}) where {N}
    p = reshape(parent(A), shp)
    return AxisArray(p, reshape_axes(naxes(A, Val(N)), axes(p)))
end

for f in (:sort, :sort!)
    @eval function Base.$f(A::AxisArray; dims, kwargs...)
        p = Base.$f(parent(A); dims=dims, kwargs...)
        return AxisArray(p, map(assign_indices, axes(A), axes(p)))
    end

    # Vector case
    @eval function Base.$f(A::AxisArray{T,1}; kwargs...) where {T}
        p = Base.$f(parent(A); kwargs...)
        return AxisArray(p, map(assign_indices, axes(A), axes(p)))
    end
end

###
### init_array
###

_length(x::Integer) = x
_length(x) = length(x)

#= TODO replace this with what's in VectorizationBase.jl
function init_array(::Type{T}, init::ArrayInitializer, axs::NTuple{N,Any}) where {T,N}
    create_static_array = true
    for i in 1:N
        is_static(getfield(axs, i)) || return Array{T,N}(init, map(_length, axs))
    end
    return MArray{Tuple{map(_length, axs)...},T,N}(init)
end
=#
function init_array(::Type{T}, init::ArrayInitializer, axs::NTuple{N,Any}) where {T,N}
    return Array{T,N}(init, map(_length, axs))
end

#=
function static_init_array(::Type{T}, init::ArrayInitializer, sz::NTuple{N,Any}) where {T,N}
    return
end

function fixed_init_array(::Fixed, ::Type{T}, init::ArrayInitializer, sz::NTuple{N,Any}) where {T,N}
    return
end

# TODO
# Currently, the only dynamic array we support is Vector, eventually it would be
# nice if we could support >1 dimensions being dynamic
function init_array(::Dynamic, ::Type{T}, init::ArrayInitializer, sz::NTuple{N,Any}) where {T,N}
    return Array{T,N}(undef, map(_length, sz))
end
=#

Base.dataids(A::AxisArray) = Base.dataids(parent(A))

function Base.zeros(::Type{T}, axs::Tuple{Vararg{<:AbstractAxis}}) where {T}
    return initialize_axis_array(zeros(T, map(length, axs)), axs)
end

function Base.falses(axs::Tuple{Vararg{<:AbstractAxis}})
    return initialize_axis_array(falses(map(length, axs)), axs)
end

function Base.fill(x, axs::Tuple{Vararg{<:AbstractAxis}})
    return initialize_axis_array(fill(x, map(length, axs)), axs)
end

function Base.reshape(A::AbstractArray, shp::Tuple{<:AbstractAxis,Vararg{<:AbstractAxis}})
    p = reshape(parent(A), map(length, shp))
    axs = reshape_axes(naxes(shp, Val(length(shp))), axes(p))
    return initialize_axis_array(p, axs)
end

# FIXME
# When I use Val(N) on the tuple the it spits out many lines of extra code.
# But without it it loses inferrence
function Base.reinterpret(::Type{Tnew}, A::AxisArray{Told,N}) where {Tnew,Told,N}
    p = reinterpret(Tnew, parent(A))
    axs = ntuple(N) do i
        resize_last(axes(A, i), size(p, i))
    end
    return AxisArray(p, axs)
end

function Base.reverse(x::AxisArray{T,N}; dims::Integer) where {T,N}
    p = reverse(parent(x), dims=dims)
    axs = ntuple(Val(N)) do i
        if i in dims
            reverse_keys(axes(x, i), axes(p, i))
        else
            assign_indices(axes(x, i), axes(p, i))
        end
    end
    return AxisArray(p, axs)
end

Base.has_offset_axes(A::AxisArray) = Base.has_offset_axes(parent(A))

###
### Indexing
###
function unsafe_view(A, inds::Tuple{Vararg{<:Integer}})
    return @inbounds(Base.view(parent(A), apply_offsets(A, inds)...))
end

function unsafe_view(A, inds::Tuple)
    p = @inbounds(Base.view(parent(A), apply_offsets(A, inds)...))
    return unsafe_reconstruct(A, p; axes=to_axes(A, inds))
end

function unsafe_dotview(A, inds::Tuple{Vararg{<:Integer}})
    return @inbounds(Base.dotview(parent(A), _sub_offset(A, inds)...))
end

function unsafe_dotview(A, inds::Tuple)
    p = @inbounds(Base.dotview(parent(A), _sub_offset(A, inds)...))
    return AxisArray(p, to_axes(A, axes(p)))
end

@propagate_inbounds function Base.getindex(A::AxisArray, args...)
    return ArrayInterface.getindex(A, args...)
end

Base.getindex(A::AxisArray, ::Ellipsis) = A

@propagate_inbounds function Base.setindex!(A::AxisArray, val, args...)
    return ArrayInterface.setindex!(A, val, args...)
end

for (unsafe_f, f) in (
    (:unsafe_view, :view),
    (:unsafe_dotview, :dotview))
    @eval begin
        @propagate_inbounds function Base.$f(A::AxisArray, args...)
            return unsafe_view(A, to_indices(A, args))
        end
    end
end

###
### Math
###
for f in (:sum!, :prod!, :maximum!, :minimum!)
    for (A,B) in (
        (AxisArray, AbstractArray),
        (AbstractArray, AxisArray),
        (AxisArray, AxisArray))
        @eval begin
            function Base.$f(a::$A, b::$B)
                Base.$f(parent(a), parent(b))
                return a
            end
        end
    end
end

for f in (:cumsum, :cumprod)
    @eval function Base.$f(a::AxisArray; dims, kwargs...)
        p = Base.$f(parent(a); dims=dims, kwargs...)
        return AxisArray(p, map(assign_indices, axes(a), axes(p)))
    end

    # Vector case
    @eval function Base.$f(a::AxisArray{T,1}; kwargs...) where {T}
        p = Base.$f(parent(a); kwargs...)
        return AxisArray(p, map(assign_indices, axes(a), axes(p)))
    end
end

function Base.unsafe_convert(::Type{Ptr{T}}, x::AxisArray{T}) where {T}
    return Base.unsafe_convert(Ptr{T}, parent(x))
end

function Base.read!(io::IO, a::AxisArray)
    read!(io, parent(a))
    return a
end

Base.write(io::IO, a::AxisArray) = write(io, parent(a))

function Base.empty!(a::AxisArray)
    for axis in axes(a)
        if !can_set_length(axis)
            error("Cannot perform `empty!` on AxisArray that has an axis with a fixed size.")
        end
    end

    for axis in axes(a)
        empty!(axis)
    end
    empty!(parent(a))
    return a
end

function Base.convert(::Type{T}, A::AbstractArray) where {T<:AxisArray}
    if A isa T
        return A
    else
        return T(A)
    end
end

const ReinterpretAxisArray{T,N,S,A<:AxisArray{S,N}} = ReinterpretArray{T,N,S,A}

function Base.axes(A::ReinterpretAxisArray{T,N,S}) where {T,N,S}
    paxs = axes(parent(A))
    axis_1 = first(paxs)
    len = div(length(axis_1) * sizeof(S), sizeof(T))
    return tuple(resize_last(axis_1, len), tail(paxs)...)
end

function Base.collect(A::AxisArray{T,N}) where {T,N}
    p = similar(parent(A), size(A))
    copyto!(p, A)
    axs = map(unsafe_reconstruct,  axes(A), axes(p))
    return initialize_axis_array(p, axs)
end

#=
function size(a::ReinterpretArray{T,N,S} where {N}) where {T,S}
    psize = size(a.parent)
    size1 = div(psize[1]*sizeof(S), sizeof(T))
    tuple(size1, tail(psize)...)
end
=#

@inline function Base.selectdim(A::AxisArray{T,N}, d::Integer, i) where {T,N}
    axs = ntuple(N) do dim_i
        if dim_i == d
            i
        else
            (:)
        end
    end
    return view(A, axs...)
end

"""
    diag(M::AxisMatrix, k::Integer=0; dim::Val=Val(1))

The `k`th diagonal of an `AxisMatrixMatrix`, `M`. The keyword argument
`dim` specifies which which dimension's axis to preserve, with the default being
the first dimension. This can be change by specifying `dim=Val(2)` instead.

```jldoctest
julia> using AxisIndices, LinearAlgebra

julia> A = AxisArray([1 2 3; 4 5 6; 7 8 9], ["a", "b", "c"], [:one, :two, :three]);

julia> keys.(axes(diag(A)))
(["a", "b", "c"],)

julia> keys.(axes(diag(A, 1; dim=Val(2))))
([:one, :two],)

```
"""
function LinearAlgebra.diag(M::AxisArray, k::Integer=0; dim::Val{D}=Val(1)) where {D}
    p = diag(parent(M), k)
    return AxisArray(p, (StaticRanges.shrink_last(axes(M, D), axes(p, 1)),))
end

"""
    inv(M::AxisMatrix)

Computes the inverse of an `AxisMatrixMatrix`
```jldoctest
julia> using AxisIndices, LinearAlgebra

julia> M = AxisArray([2 5; 1 3], ["a", "b"], [:one, :two]);

julia> keys.(axes(inv(M)))
([:one, :two], ["a", "b"])

```
"""
function Base.inv(A::AxisArray)
    p = inv(parent(A))
    axs = (assign_indices(axes(A, 2), axes(p, 1)), assign_indices(axes(A, 1), axes(A, 2)))
    return AxisArray(p, axs)
end

for f in (
    :(Base.transpose),
    :(Base.adjoint),
    :(LinearAlgebra.pinv))
    @eval begin
        function $f(A::AxisArray)
            p = $f(parent(A))
            return AxisArray(p, permute_axes(A, p))
        end
    end
end

#=

axs = Base.Iterators.ProductIterator{Tuple{Base.OneTo{Int64},Tuple{Colon}}}((Base.OneTo(4), (Colon(),)))

axs = Base.Iterators.ProductIterator{Tuple{Base.OneTo{Int64},Tuple{Colon}}}((Base.OneTo(4), (Colon(),)))
collect(Base.Iterators.ProductIterator(axes(A)))

function Base.sortslices(A::AxisArray; dims, kwargs...)
    itspace = Base.compute_itspace(parent(A), ifelse(dims isa Val, dims, Val(dims)))
    vecs = map(its->view(parent(A), its...), itspace)
    p = sortperm(vecs; kwargs...)
    B = similar(A)
    for (x, its) in zip(p, itspace)
        B[its...] = vecs[x]
    end
    return B
end

function compute_itspace(A, ::Val{dims}) where {dims}
    N = ndims(A)
    negdims = filter(i-> !(i in dims), 1:N)
    axs = Iterators.product(ntuple(DimSelector{dims}(A), N)...)
    vec(permutedims(collect(axs), (dims..., negdims...)))
end
=#

for f in (:map, :map!)
    # Here f::F where {F} is needed to avoid ambiguities in Julia 1.0
    @eval begin
        function Base.$f(f::F, a::AbstractArray, b::AxisArray, cs::AbstractArray...) where {F}
            return AxisArray(
                $f(f, parent(a), parent(b), parent.(cs)...),
                Broadcast.combine_axes(a, b, cs...,)
            )
        end

        function Base.$f(f::F, a::AxisArray, b::AxisArray, cs::AbstractArray...) where {F}
            return AxisArray(
                $f(f, parent(a), parent(b), parent.(cs)...),
                Broadcast.combine_axes(a, b, cs...,)
            )
        end

        function Base.$f(f::F, a::AxisArray, b::AbstractArray, cs::AbstractArray...) where {F}
            return AxisArray(
                $f(f, parent(a), parent(b), parent.(cs)...),
                Broadcast.combine_axes(a, b, cs...,)
            )
        end
    end
end

Base.map(f, A::AxisArray) = AxisArray(map(f, parent(A)), axes(A))

# We can't just make a type alias for mapped array types because this would require
# multiple calls to combine_axes for multi-mapped types for every axes call. It also
# would require overloading a bunch of other methods to ensure they work correctly
# (e.g., getindex, setindex!, view, show, etc...)
#
# We can't directly overload the head of each method because data::AbstractArray....
# is too similar to Union{AxisArray,AbstractArray} so we only specialize
# on method heads that handle all AxisArray subtypes. Therefore, including
# any other array type will miss these specific methods.

function MappedArrays.mappedarray(f, data::AxisArray)
    return AxisArray(mappedarray(f, parent(data)), axes(data))
end

function MappedArrays.mappedarray(::Type{T}, data::AxisArray) where T
    return AxisArray(mappedarray(T, parent(data)), axes(data))
end

function MappedArrays.mappedarray(f, data::AxisArray...)
    return AxisArray(
        mappedarray(f, map(unwrap_broadcasted, data)...),
        Broadcast.combine_axes(data...)
    )
end

function MappedArrays.mappedarray(::Type{T}, data::AxisArray...) where T
    return AxisArray(
        mappedarray(T, map(unwrap_broadcasted, data)...),
        Broadcast.combine_axes(data...)
    )
end

# These needed to have the additional ::Function defined to avoid ambiguities
function MappedArrays.mappedarray(f, finv::Function, data::AxisArray)
    return AxisArray(mappedarray(f, finv, parent(data)), axes(data))
end

function MappedArrays.mappedarray(f, finv::Function, data::AxisArray...)
    return AxisArray(
        mappedarray(f, finv, map(unwrap_broadcasted, data)...),
        Broadcast.combine_axes(data...)
    )
end

function MappedArrays.mappedarray(::Type{T}, finv::Function, data::AxisArray...) where T
    return AxisArray(
        mappedarray(T, finv, map(unwrap_broadcasted, data)...),
        Broadcast.combine_axes(data...)
    )
end

function MappedArrays.mappedarray(f, ::Type{Finv}, data::AxisArray...) where Finv
    return AxisArray(
        mappedarray(f, Finv, map(unwrap_broadcasted, data)...),
        Broadcast.combine_axes(data...)
    )
end

function MappedArrays.mappedarray(::Type{T}, ::Type{Finv}, data::AxisArray...) where {T,Finv}
    return AxisArray(
        mappedarray(T, Finv, map(unwrap_broadcasted, data)...),
        Broadcast.combine_axes(data...)
    )
end

#
#    AxisArrayStyle{S}
#
# This is a `BroadcastStyle` for AxisArray's It preserves the dimension
# names. `S` should be the `BroadcastStyle` of the wrapped type.
struct AxisArrayStyle{S <: BroadcastStyle} <: AbstractArrayStyle{Any} end
AxisArrayStyle(::S) where {S} = AxisArrayStyle{S}()
AxisArrayStyle(::S, ::Val{N}) where {S,N} = AxisArrayStyle(S(Val(N)))
AxisArrayStyle(::Val{N}) where N = AxisArrayStyle{DefaultArrayStyle{N}}()
function AxisArrayStyle(a::BroadcastStyle, b::BroadcastStyle)
    inner_style = BroadcastStyle(a, b)

    # if the inner_style is Unknown then so is the outer-style
    if inner_style isa Unknown
        return Unknown()
    else
        return AxisArrayStyle(inner_style)
    end
end

function Base.BroadcastStyle(::Type{T}) where {T<:AxisArray}
    return AxisArrayStyle{typeof(BroadcastStyle(parent_type(T)))}()
end

Base.BroadcastStyle(::AxisArrayStyle{A}, ::AxisArrayStyle{B}) where {A, B} = AxisArrayStyle(A(), B())
#Base.BroadcastStyle(::AxisArrayStyle{A}, b::B) where {A, B} = AxisArrayStyle(A(), b)
#Base.BroadcastStyle(a::A, ::AxisArrayStyle{B}) where {A, B} = AxisArrayStyle(a, B())
Base.BroadcastStyle(::AxisArrayStyle{A}, b::DefaultArrayStyle) where {A} = AxisArrayStyle(A(), b)
Base.BroadcastStyle(a::AbstractArrayStyle{M}, ::AxisArrayStyle{B}) where {B,M} = AxisArrayStyle(a, B())

#
#    unwrap_broadcasted
#
# Recursively unwraps `AxisArray`s and `AxisArrayStyle`s.
# replacing the `AxisArray`s with the wrapped array,
# and `AxisArrayStyle` with the wrapped `BroadcastStyle`.
function unwrap_broadcasted(bc::Broadcasted{AxisArrayStyle{S}}) where S
    return Broadcasted{S}(bc.f, map(unwrap_broadcasted, bc.args))
end
unwrap_broadcasted(a::AxisArray) = parent(a)
unwrap_broadcasted(x) = x

# We need to implement copy because if the wrapper array type does not support setindex
# then the `similar` based default method will not work
function Broadcast.copy(bc::Broadcasted{AxisArrayStyle{S}}) where S
    return AxisArray(copy(unwrap_broadcasted(bc)), Broadcast.combine_axes(bc.args...))
end

for f in (:zero, :one)
    @eval begin
        function Base.$f(a::AxisArray)
            p = Base.$f(parent(a))
            return AxisArray(p, map(assign_indices, axes(a), axes(p)))
        end
    end
end

drop_axes(x::AbstractArray, d::Int) = drop_axes(x, (d,))
drop_axes(x::AbstractArray, d::Tuple) = drop_axes(x, dim(dimnames(x), d))
drop_axes(x::AbstractArray, d::Tuple{Vararg{Int}}) = drop_axes(axes(x), d)
drop_axes(x::Tuple{Vararg{<:Any}}, d::Int) = drop_axes(x, (d,))
drop_axes(x::Tuple{Vararg{<:Any}}, d::Tuple) = _drop_axes(x, d)
_drop_axes(x, y) = select_axes(x, dropinds(x, y))

select_axes(x::AbstractArray, d::Tuple) = select_axes(x, dims(dimnames(x), d))
select_axes(x::AbstractArray, d::Tuple{Vararg{Int}}) = map(i -> axes(x, i), d)
select_axes(x::Tuple, d::Tuple) = map(i -> getfield(x, i), d)

dropinds(x, y) = _dropinds(x, y)
Base.@pure @inline function _dropinds(x::Tuple{Vararg{Any,N}}, dims::NTuple{M,Int}) where {N,M}
    out = ()
    for i in 1:N
        cnd = true
        for j in dims
            if i === j
                cnd = false
                break
            end
        end
        if cnd
            out = (out..., i)
        end
    end
    return out::NTuple{N - M, Int}
end

function Base.dropdims(a::AxisArray; dims)
    return AxisArray(dropdims(parent(a); dims=dims), drop_axes(a, dims))
end

for (X,Y) in (
    (:(Base.Indices),:(Tuple{Vararg{<:AbstractAxis}})),
    (:(Tuple{Vararg{<:AbstractAxis}}),:(Base.Indices)),
    (:(Tuple{Vararg{<:AbstractAxis}}),:(Tuple{Vararg{<:AbstractAxis}})))
    @eval begin
        function Base.promote_shape(a::$X, b::$Y)
            if length(a) < length(b)
                return Base.promote_shape(b, a)
            end
            for i=1:length(b)
                if length(a[i]) != length(b[i])
                    throw(DimensionMismatch("dimensions must match: a has dims $a, b has dims $b, mismatch at $i"))
                end
            end
            for i=length(b)+1:length(a)
                if length(a[i]) != 1
                    throw(DimensionMismatch("dimensions must match: a has dims $a, must have singleton at dim $i"))
                end
            end
            return a
        end
    end
end


text:
#Code written by Skipper af Rosenborg and Esben Bager
using StatsBase
using DataFrames
using CSV
using JuMP
include("SupportFunction.jl")
include("DataLoad.jl")
println("Leeeeroooy Jenkins")

#Aggregate data for each of i gamma values
for i = 1:16
    #set gamma value
    gammaArr = linspace(0,4,16)

    #Set name and path of dataset
    dataSet = "VIXBestRegulation740-790/"*string(i)*"/"
    path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/CEO/RFR/"*dataSet
    cd(path)

    #initilise arrays
    totalMatrix = zeros(844,12)
    forecastMatrix = zeros(844,22)

    for i = 0:84
        #Load data
        tempMatrix = CSV.read(path*string(i)*"_returnPvalueOutcome1to200.csv",nullable=false, header=false, datarow=1, types=fill(Float64,4))
        tempForecast = CSV.read(path*string(i)*"_forecastsAndErrors1to200.csv",nullable=false, header=false, datarow=1, types=fill(Float64,22))
        if i == 0
            totalMatrix = Array(tempMatrix)
            forecastMatrix = Array(tempForecast)
        else
            totalMatrix[6:end,2:end] += Array(tempMatrix[6:end,2:end])
            forecastMatrix[6:end,1:end] += Array(tempForecast[6:end,1:end])
        end
    end

    writedlm(path*"VIXTimeTABestRegForecast"*string(round(gammaArr[i],3))*".csv", forecastMatrix, ",")
    writedlm(path*"VIXTimeTABestReg"*string(round(gammaArr[i],3))*".csv",totalMatrix,",")

    weightMatrix = zeros(846,11)
    for i = 1:1
        weightMatrix = zeros(846,11)
        for j = 0:84
            tempMatrix = CSV.read(path*"wightsCEO_Model"*string(i)*"_"*string(j)*".csv",nullable=false, header=false, datarow=1, types=fill(Float64,11))
            if j == 0
                weightMatrix = Array(tempMatrix)
            else
                weightMatrix += Array(tempMatrix)
            end
        end
        writedlm(path*"weightsModel"*string(round(gammaArr[i],3))*".csv",weightMatrix[1:844,:],",")
    end
end

#Aggregate return data from i gamma periods into 1 file
totalArray = zeros(844,34)
tempString = "Perfect,1/N"
for i = 1:16
    gammaArr = linspace(0,4,16)
    dataSet = "VIXBestRegulation740-790/"*string(i)*"/"
    path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/CEO/RFR/"*dataSet

    tempArr = CSV.read(path*"VIXTimeTABestReg"*string(round(gammaArr[i],3))*".csv", types = fill(Float64,4),
        datarow = 1, header=false, nullable=false)

    if i == 1
        totalArray[:,1:2] = Array{Float64}(tempArr)[:,1:2]
    end

    totalArray[:, 1+i*2:2+i*2] = Array{Float64}(tempArr)[:,3:4]

    tempString = tempString*",Return "*string(round(gammaArr[i],3))*",Pvalue "*string(round(gammaArr[i],3))
end
path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/CEO/RFR/"*"VIXBestRegulation740-790/"

f = open(path*"TotalReturnArray.csv", "w")
write(f, tempString*"\n")
writecsv(f, totalArray)
close(f)

max_stars_repo_name: skipperafrosenborg/Thesis

max_issues_repo_path: ML/Code to hand in/DataAnalysis_CEO.jl

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_path: ML/Code to hand in/DataAnalysis_CEO.jl

max_forks_repo_name: skipperafrosenborg/Thesis


text:
module TDD
export Interval
using DataStructures
# Write your package code here.
include("graph.jl")
include("interval.jl")
end


text:
# This file was generated, do not modify it. # hide
linear = (@load LinearRegressor pkg=MLJLinearModels)()
knn = (@load KNNRegressor)()

tree_booster = (@load EvoTreeRegressor)()
forest = (@load RandomForestRegressor pkg=DecisionTree)()
svm = (@load SVMRegressor)()

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
"""
Module that can create QR codes as data or images using `qrcode` or `exportqrcode`.
"""
module QRCode

export Mode, Numeric, Alphanumeric, Byte
export ErrCorrLevel, Low, Medium, Quartile, High
export getmode, getversion, qrcode, exportqrcode

using ImageCore
using FileIO

"""
Abstract type that groups the three supported encoding modes `Numeric`,
`Alphanumeric` and `Byte`.
"""
abstract type Mode end
"""
Encoding mode for messages composed of digits only.
"""
struct Numeric <: Mode end
"""
Encoding mode for messages composed of digits, characters `A`-`Z` (capital only)
, space and `%` `*` `+` `-` `.` `/` `:` `\$`.
"""
struct Alphanumeric <: Mode end
"""
Encoding mode for messages composed of ISO 8859-1 or UTF-8 characters.
"""
struct Byte <: Mode end
# struct Kanji <: Mode end

"""
Abstract type that groups the four error correction levels `Low`, `Medium`,
`Quartile` and `High`.
"""
abstract type ErrCorrLevel end
"""
Error correction level that can restore up to 7% of missing codewords.
"""
struct Low <: ErrCorrLevel end
"""
Error correction level that can restore up to 15% of missing codewords.
"""
struct Medium <: ErrCorrLevel end
"""
Error correction level that can restore up to 25% of missing codewords.
"""
struct Quartile <: ErrCorrLevel end
"""
Error correction level that can restore up to 30% of missing codewords.
"""
struct High <: ErrCorrLevel end

include("tables.jl")
include("errorcorrection.jl")
include("matrix.jl")

using .Polynomial

"""
    getmode(message::AbstractString)

Return the encoding mode of `message`, either `Numeric()`, `Alphanumeric()`
or `Byte()`.

# Examples
```jldoctest
julia> getmode("HELLO WORLD")
Alphanumeric()
```
"""
function getmode(message::AbstractString)::Mode
    if all([isdigit(c) for c in message])
        return Numeric()
    elseif all([haskey(alphanumeric, c) for c in message])
        return Alphanumeric()
    else
        return Byte()
    end
end

"""
    getversion(message::AbstractString, mode::Mode, level::ErrCorrLevel)

Return the version of the QR code, between 1 and 40.

```jldoctest
julia> getversion("Hello World!", Alphanumeric(), High())
2
```
"""
function getversion(message::AbstractString, mode::Mode, level::ErrCorrLevel)
    cc = characterscapacity[(level, mode)]
    return findfirst(v->v >= lastindex(message), cc)
end

"""
    getcharactercountindicator(msglength::Int, version::Int, mode::Mode)

Return the bits for character count indicator.
"""
function getcharactercountindicator(msglength::Int,
                                    version::Int,
                                    mode::Mode)::BitArray{1}
    if 1 <= version <= 9
        i = 1
    elseif version <= 26
        i = 2
    else
        i = 3
    end
    cclength = charactercountlength[mode][i]
    indicator = BitArray(reverse(digits(msglength, base = 2, pad = cclength)))

    return indicator
end

"""
    encodedata(message::AbstractString, ::Mode)

Encode the message with the given mode.
"""
function encodedata(message::AbstractString, ::Numeric)::BitArray{1}
    l = length(message)
    chunks = [SubString(message, i, min(i + 2, l)) for i in 1:3:l]

    function toBin(chunk::SubString)::BitArray
        pad = 1 + 3 * length(chunk)
        n = parse(Int, chunk)
        return  BitArray(reverse(digits(n, base = 2, pad = pad)))
    end

    binchunks = map(toBin, chunks)
    return vcat(binchunks...)
end

function encodedata(message::AbstractString, ::Alphanumeric)::BitArray{1}
    l = length(message)
    chunks = [SubString(message, i, min(i + 1, l)) for i in 1:2:l]

    function toBin(s::SubString)::BitArray{1}
        if length(s) == 1
            n = alphanumeric[s[1]]
            return  BitArray{1}(reverse(digits(n, base = 2, pad = 6)))
        else
            n = 45 * alphanumeric[s[1]] + alphanumeric[s[2]]
            return  BitArray{1}(reverse(digits(n, base = 2, pad = 11)))
        end
    end

    binchunks = map(toBin, chunks)
    return vcat(binchunks...)
end

function encodedata(message::AbstractString, ::Byte)::BitArray{1}
    bytes = Array{UInt8}(message)
    bin = map(int2bitarray, Array{Int,1}(bytes))
    return vcat(bin...)
end

"""
    int2bitarray(n::Int)

Encode an integer into a `BitArray`.
"""
function int2bitarray(n::Int)::BitArray{1}
    return  BitArray(reverse(digits(n, base = 2, pad = 8)))
end

"""
    padencodedmessage(data::BitArray{1}, requiredlentgh::Int)

Pad the encoded message.
"""
function padencodedmessage(data::BitArray{1}, requiredlentgh::Int)
    # Add up to 4 zeros to terminate the message
    data = vcat(data, falses(min(4, requiredlentgh - length(data))))

    # Add zeros to make the length a multiple of 8
    if length(data) % 8 != 0
        data = vcat(data, falses(8 - length(data) % 8))
    end

    # Add the repeated pattern until reaching required length
    pattern = BitArray{1}([1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1])
    pad = repeat(pattern, ceil(Int, requiredlentgh - length(data) / 8))
    data = vcat(data, pad[1:requiredlentgh - length(data)])

    return data
end

"""
    makeblocks(bits::BitArray{1}, nb1::Int, dc1::Int, nb2::Int, dc2::Int)

Divide the encoded message into 1 or 2 blocks with `nbX` blocks in group `X` and
`dcX` codewords per block.
"""
function makeblocks(bits::BitArray{1},
                    nb1::Int,
                    dc1::Int,
                    nb2::Int,
                    dc2::Int)::Array{BitArray{2},1}
    array = reshape(bits, (8, length(bits) ÷ 8))
    i = 1
    blocks = BitArray{2}[]
    for n in 1:nb1
        push!(blocks, array[:, i:i + dc1 - 1])
        i += dc1
    end
    for n in 1:nb2
        push!(blocks, array[:, i:i + dc2 - 1])
        i += dc2
    end
    return blocks
end

"""
    geterrcorrblock(block::BitArray{2}, ncodewords::Int)

Return the error correction blocks, with `ncodewords` codewords per block.
"""
function geterrcorrblock(block::BitArray{2}, ncodewords::Int)::BitArray{2}
    # Helper functions
    bitarray2ints(b) = reduce((acc, n)->2 * acc + n, b, init = 0, dims = 1)
    array2poly(a) =  Poly(reverse(a[:]))

    poly = array2poly(bitarray2ints(block))
    ecpoly = geterrorcorrection(poly, ncodewords)
    ecarray = map(int2bitarray, reverse(ecpoly.coeff))
    ecbits = foldl(vcat, ecarray, init = BitArray{1}())
    return reshape(ecbits, (8, length(ecbits) ÷ 8))
end

"""
    interleave(blocks::Array{BitArray{2},1}, ecblocks::Array{BitArray{2},1},
               ncodewords::Int, nb1::Int, dc1::Int, nb2::Int, dc2::Int,
               version::Int)

Mix the encoded data blocks and error correction blocks.
"""
function interleave( blocks::Array{BitArray{2},1}
                   , ecblocks::Array{BitArray{2},1}
                   , ncodewords::Int, nb1::Int
                   , dc1::Int
                   , nb2::Int
                   , dc2::Int
                   , version::Int
                   )::BitArray{1}

    data = BitArray{1}()

    slice(i) = (acc, block) -> vcat(acc, block[:, i])

    # Encoded data
    for i in 1:dc1
        data = foldl(slice(i), blocks, init = data)
    end
    if dc2 > dc1
        data = foldl(slice(dc2), blocks[nb1 + 1 : nb1 + nb2], init = data)
    end

    # Error correction data
    for i in 1:ncodewords
        data = foldl(slice(i), ecblocks, init = data)
    end

    # Extra padding
    data = vcat(data, falses(remainderbits[version]))

    return data
end

"""
    qrcode(message::AbstractString, eclevel = Medium(); compact = false)

Create a `BitArray{2}` with the encoded `message`, with `true` (`1`) for the black
areas and `false` (`0`) as the white ones. If `compact` is `false`, white space
is added around the QR code.

The error correction level `eclevel` can be picked from four values: `Low()`
(7% of missing data can be restored), `Medium()` (15%), `Quartile()` (25%) or
`High()` (30%). Higher levels make denser QR codes.
"""
function qrcode( message::AbstractString
               , eclevel::ErrCorrLevel = Medium()
               ; compact::Bool = false )
    # Determining QR code mode and version
    mode = getmode(message)
    version = getversion(message, mode, eclevel)

    # Mode indicator: part of the encoded message
    modeindicator = modeindicators[mode]

    # Character count: part of the encoded message
    ccindicator = getcharactercountindicator(lastindex(message), version, mode)

    # Encoded data: main part of the encoded message
    encodeddata = encodedata(message, mode)

    # Getting parameters for the error correction
    # Number of error correction codewords per block, number of blocks in
    # group 1/2, number of data codewords per block in group 1/2
    ncodewords, nb1, dc1, nb2, dc2 = ecblockinfo[eclevel][version, :]
    requiredbits = 8 * (nb1 * dc1 + nb2 * dc2)

    # Pad encoded message before error correction
    encoded = vcat(modeindicator, ccindicator, encodeddata)
    encoded = padencodedmessage(encoded, requiredbits)

    # Getting error correction codes
    blocks = makeblocks(encoded, nb1, dc1, nb2, dc2)
    ecblocks = map(b->geterrcorrblock(b, ncodewords), blocks)

    # Interleave code blocks
    data = interleave(blocks, ecblocks, ncodewords, nb1, dc1, nb2, dc2, version)

    # Generate qr code matrix, masks and fill it
    matrix = emptymatrix(version)
    masks = makemasks(matrix)
    matrix = placedata!(matrix, data)

    # Pick the best mask
    candidates = map(enumerate(masks)) do (i, m)
        i - 1, xor.(matrix, m)
    end
    mask, matrix = first(sort(candidates, by = penalty ∘ last))

    # Format and version information
    matrix = addformat(matrix, mask, version, eclevel)

    if compact
        return matrix
    else
        background = falses(size(matrix) .+ (8, 8))
        background[5:end-4, 5:end-4] = matrix
        return background
    end
end

"""
    exportqrcode( message::AbstractString
                , path = "qrcode.png"
                , eclevel = Medium()
                ; targetsize = 5
                , compact = false )

Create a `PNG` file with the encoded `message` of approximate size `targetsize`
cm. If `compact` is `false`, white space is added around the QR code.

The error correction level `eclevel` can be picked from four values: `Low()`
(7% of missing data can be restored), `Medium()` (15%), `Quartile()` (25%) or
`High()` (30%). Higher levels make denser QR codes.
"""
function exportqrcode( message::AbstractString
                     , path::AbstractString = "qrcode.png"
                     , eclevel::ErrCorrLevel = Medium()
                     ; targetsize::Int = 5
                     , compact::Bool = false )

    matrix = qrcode(message, eclevel, compact = compact)

    if !endswith(path, ".png")
        path = "$path.png"
    end

    # Seems the default setting is 72 DPI
    pixels = size(matrix, 1)
    scale = ceil(Int, 72 * targetsize / 2.45 / pixels)
    matrix = kron(matrix, trues((scale, scale)))

    save(path, BitArray(.! matrix))
end

end # module


text:
using JudiLing
using Test
using CSV
using DataFrames

@testset "make cue matrix for latin" begin
    latin = DataFrame(CSV.File(joinpath("data", "latin_mini.csv")))
    cue_obj = JudiLing.make_cue_matrix(
        latin,
        grams = 3,
        target_col = :Word,
        tokenized = false,
        keep_sep = false,
    )
    JudiLing.make_Yt_matrix(2, latin, cue_obj.f2i, target_col = :Word)
end


text:
@testset "recession_bands" begin
    p1 = Plot(scatter(x=1:10, y=rand(10)))
    @test PlotlyBase._recession_band_shapes(p1) === nothing

    _path = joinpath(@__DIR__, "data", "us_manu_unemp.csv")
    _with_dates_raw = readdlm(_path, ',')
    with_dates = _with_dates_raw[2:end, :]
    with_dates[:, 1] = map(Date, with_dates[:, 1])

    p2 = Plot(scatter(x=with_dates[:, 1], y=with_dates[:, 2]))
    @test PlotlyBase._recession_band_shapes(p2) !== nothing
    @test length(PlotlyBase._recession_band_shapes(p2)) == 7

    p3 = Plot(scatter(x=with_dates[:, 1], y=with_dates[:, 3]))
    @test PlotlyBase._recession_band_shapes(p3) !== nothing
    @test length(PlotlyBase._recession_band_shapes(p3)) == 7

    p12 = [p1; p2]
    @test PlotlyBase._recession_band_shapes(p12) !== nothing
    @test length(PlotlyBase._recession_band_shapes(p12)) == 7

    p123 = [p1; p2; p3]
    @test PlotlyBase._recession_band_shapes(p123) !== nothing
    @test length(PlotlyBase._recession_band_shapes(p123)) == 14

    add_recession_bands!(p1)
    @test length(p1.layout["shapes"]) == 0

    add_recession_bands!(p123)
    @test length(p123.layout["shapes"]) == 14

    p23 = [p2; p3]
    china_WTO_date = Date(2001, 12, 11)
    relayout!(
        p23,
        shapes=vline([china_WTO_date], line_color="red", opacity=0.4),
    )
    add_recession_bands!(p23)
    @test length(p23.layout["shapes"]) == 15

end


text:
# --------------- #
# Chebyshev Basis #
# --------------- #

function Basis(::Cheb, n::Int, a::Real, b::Real)  # chebdef
    n <= 0 && error("n must be positive")
    a >= b && error("left endpoint (a) must be less than right end point (b)")
    Basis(Cheb(), n, a, b, ChebParams(n, a, b))
end

# define methods for ChepParams type
Basis(p::ChebParams) = Basis(Cheb(), p.n, p.a, p.b)

# chebnode.m -- DONE
function nodes(p::ChebParams, nodetype=0)
    n, a, b = p.n, p.a, p.b
    s = (b-a) / 2  # 21
    m = (b+a) / 2  # 22

    if nodetype < 2  # usual nodes
        k = π*(0.5:(maximum(n) - 0.5))'  # 25
        x = m - cos(k[1:n] / n) * s  # 26
        if nodetype == 1  # Extend nodes to endpoints
            aa = x[1]  # 28
            bb = x[end]  # 29
            x = (bb*a - aa*b)/(bb-aa) + (b-a)/(bb-aa)*x  # 30
        end
    else # Lobatto nodes
        k = pi*(0:(maximum(n)-1))  # 33
        x = m - cos(k[1:n] / (n-1)) * s  # 34
    end
    x
end

function derivative_op(p::ChebParams, order=1)
    n, a, b = p.n, p.a, p.b
    if order > 0
        # TODO: figure out some caching mechanism that avoids globals
        D = Array(SparseMatrixCSC{Float64,Int64}, max(2, order))  # 49
        i = repeat(collect(1:n), outer=[1, n]); j = i'  # 50

        # 51
        inds = find((rem(i + j, 2) .== 1) & (j.>i))
        r, c = similar(inds), similar(inds)
        for i in 1:length(inds)
            r[i], c[i] = ind2sub((n, n), inds[i])
        end

        d = sparse(r, c, (4/(b-a)) * (vec(j[1, c])-1), n-1, n)  # 52
        d[1, :] ./= 2  # 53
        D[1] = d  # 54
        for ii=2:max(2, order)
            D[ii] = d[1:n-ii, 1:n-ii+1] * D[ii-1]  # 56
        end
    else
        D = Array(SparseMatrixCSC{Float64,Int64}, abs(order))  # 64
        nn = n - order  # 65
        i = (0.25 * (b - a)) ./(1:nn)  # 66
        d = sparse(vcat(1:nn, 1:nn-2), vcat(1:nn, 3:nn), vcat(i, -i[1:nn-2]),
                   nn, nn)  # 67
        d[1, 1] *= 2  # 68
        d0 = ((-1).^(0:nn-1)') .* sum(d, 1)  # 69
        D[1] = sparse(vcat(d0[1:n]', d[1:n, 1:n]))  # 70
        for ii=-2:-1:order
            ind = 1:n-ii-1
            D[-ii] = sparse(vcat(d0[ind]', d[ind, ind]) * D[-ii-1])
        end
    end
    D, ChebParams(n-order, a, b)
end

function evalbase(p::ChebParams, x=nodes(p, 1), order::Int=0, nodetype::Int=1)
    n, a, b = p.n, p.a, p.b
    minorder = min(0, order)  # 30

    # compute 0-order basis
    if nodetype == 0
        temp = ((n-0.5):-1:0.5)''  # 41
        bas = cos((pi/n)*temp.*(0:(n-1-minorder))')  # 42
    else
        bas = evalbasex(ChebParams(n-minorder, a, b), x)  # 44
    end

    if order != 0
        D = derivative_op(p, order)[1]
        B = bas[:, 1:n-order]*D[abs(order)]
    else
        B = bas
    end

    return B
end

function evalbase(p::ChebParams, x, order::AbstractVector{Int}, nodetype::Int=1)
    n, a, b = p.n, p.a, p.b
    minorder = min(0, minimum(order))  # 30
    maxorder = maximum(order)

    # compute 0-order basis
    if nodetype == 0
        temp = ((n-0.5):-1:0.5)''  # 41
        bas = cos((pi/n)*temp.*(0:(n-1-minorder))')  # 42
    else
        bas = evalbasex(ChebParams(n-minorder, a, b), x)  # 44
    end

    B = Array(Matrix{Float64}, length(order))
    if maxorder > 0 D = derivative_op(p, maxorder)[1] end
    if minorder < 0 I = derivative_op(p, minorder)[1] end

    for ii=1:length(order)
        if order[ii] == 0
            B[ii] = bas[:, 1:n]
        elseif order[ii] > 0
            B[ii] = bas[:, 1:n-order[ii]] * D[order[ii]]
        else
            B[ii] = bas[:, 1:n-order[ii]] * I[-order[ii]]
        end
    end

    return B
end

function evalbasex(p::ChebParams, x)
    n, a, b = p.n, p.a, p.b
    z = (2/(b-a)) * (x-(a+b)/2)
    m = size(z, 1)
    bas = Array(Float64, m, n);
    bas[:, 1] = 1.0
    bas[:, 2] = z
    z = 2 * z
    @inbounds for j=3:n
        @simd for i=1:m
            bas[i, j] = z[i] .* bas[i, j-1] - bas[i, j-2]
        end
    end
    bas
end


text:
function argparse(cmd_args=ARGS)
    parse_settings = ArgParseSettings(
        commands_are_required=true,
        suppress_warnings=true,
    )

    @add_arg_table parse_settings begin
        "install"
            action = :command
            help = "Installs julia version for you."
        "create"
            action = :command
            help = "Builds the playground"
        "activate"
            action = :command
            help = "Activates the playground."
        "list"
            action = :command
            help = "Lists available julia versions and playgrounds"
        "clean"
            action = :command
            help = "Deletes any dead julia-version or playground links, in case you've deleted the original folders."
        "rm"
            help = "Deletes the specifid julia-version or playground."
            action = :command
        "exec"
            action = :command
            help = "Execute a cmd inside a playground and exit."
    end

    @add_arg_table parse_settings["install"] begin
        "download"
            action = :command
            help = "The version of julia you'd like to install from the julia website."
        "link"
            action = :command
            help = "Path to an existing julia build you'd like to use with playgrounds."
        # "build"
        #     action = :command
        #     help = "The git url to clone the julia source from"
    end

    @add_arg_table parse_settings["create"] begin
        "dir"
            help = "Where to create the playground. Defaults to the current working directory (can be changed in `~/.playground/config`."
            default = ""
        "--requirements", "-r"
            help = "A REQUIRE or DECLARE file of dependencies to install into the playground."
            default = ""
        "--name", "-n"
            help = "A global name to allow activating the playground from anywhere."
            default = ""
        "--julia-version", "-j"
            help = "The version(s) of julia available to use. If multiple versions are provided the first entry will be the one used by `julia`. By default the user/system level version is used."
            default = ""
        "--req-type", "-t"
            help = "If --requirments isn't being passed a path ending in REQUIRE or DECLARE file, please specify which type is it \"REQUIRE\" or \"DECLARE\""
            arg_type = Symbol
            default = :REQUIRE
    end

    @add_arg_table parse_settings["activate"] begin
        "dir"
            help = "The path to the playground to use. This takes priority over --name."
            default = ""
        "--name", "-n"
            help = "A global name to allow activating the playground from anywhere."
            default = ""
    end

    @add_arg_table parse_settings["exec"] begin
        "cmd"
            help = "The command you would like to run inside the playground."
            required = true
        "dir"
            help = "The path to the playground to use. This takes priority over --name."
            default = ""
        "--name", "-n"
            help = "A global name to allow activating the playground from anywhere."
            default = ""
    end

    @add_arg_table parse_settings["list"] begin
        "--show-links", "-s"
            help = "Display the source path if julia-versions or playgrounds are just symlinks."
            action = :store_true
            default = false
    end

    @add_arg_table parse_settings["clean"]

    @add_arg_table parse_settings["rm"] begin
        "name"
            help = "Deletes the playground directory with the given name and the link to it."
            default = ""
        "--dir"
            help = "Deletes the provided playground directory and the link to it."
            default = ""
    end

    @add_arg_table parse_settings["install"]["download"] begin
        "version"
            help = "The release version available to download at http://julialang.org/downloads/"
            arg_type = VersionNumber
            required = true
        "--labels", "-l"
            help = "Extra labels to apply to the new julia versions."
            arg_type = AbstractString
            nargs = '*'
    end

    @add_arg_table parse_settings["install"]["link"] begin
        "dir"
            help = "The path to a julia executable you'd like to be made available to playgrounds."
            default = ""
        "--labels", "-l"
            help = "Extra labels to apply to the new julia versions."
            arg_type = AbstractString
            nargs = '*'
    end

    # @add_arg_table parse_settings["install"]["build"] begin
    #     "url"
    #         help = "The git url to clone the julialang source from. Defaults to https://github.com/JuliaLang/julia.git"
    #         default = ""
    #     "revision"
    #         help = "The revision to checkout prior to building julia. Defaults to origin/master"
    #         default = ""
    #     "--labels", "-l"
    #         help = "Extra labels to apply to the new julia versions."
    #         arg_type = AbstractString
    #         nargs = '*'
    # end

    return parse_args(cmd_args, parse_settings)
end


text:
module Otm

using ..FEA,
      LinearAlgebra

export Otimizador,
       otimizar!,
       gerar_trelica_hiperconectada

"""
    gerar_trelica_hiperconectada(Lx::Real,
                                 Ly::Real,
                                 nx::Integer,
                                 ny::Integer,
                                 material::Material,
                                 nivel::Integer=2,
                                 area::Real=1)::Estrutura

Gera uma treliça hiperconectada de domínio estendido retangular.

# Argumentos
- `Lx::Real`: Comprimento ao longo de x.
- `Ly::Real`: Comprimento ao longo de y.
- `nx::Integer`: Número de nós na direção x.
- `ny::Integer`: Número de nós na direção y.
- `material::Material`: Material de todos os elementos.
- `nivel::Integer=2`: Nível de conectividade dos nós.
- `area::Real=1`: Área inicial das seções dos elementos.
"""
function gerar_trelica_hiperconectada(Lx::Real,
                                      Ly::Real,
                                      nx::Integer,
                                      ny::Integer,
                                      material::Material,
                                      nivel::Integer=2,
                                      area::Real=1)::Estrutura
    # Distância entre nós na direção x.
    lx = Lx/(nx - 1)
    # Distância entre nós na direção y.
    ly = Ly/(ny - 1)
    # Número total de nós na estrutura.
    n = nx * ny

    # Criação da matriz que contém as posições de cada nó.
    m = reshape(1:n, (nx, ny))'
    # Criação dos nós da estrutura.
    nos::Vector{No} = []
    k = 1
    for i=1:ny
        y = Ly - (i - 1)*ly
        for j=1:nx
            x = (j - 1)*lx
            push!(nos, No(k, x, y))
            k += 1
        end
    end

    # Criação da treliça hiperconectada sem que haja sobreposição de barras colineares.
    elementos::Vector{Elemento} = []
    c = 1
    for i=1:ny
        for j=1:nx
            no_ref = nos[m[i, j]]
            @info "Gerando treliça: $(m[i, j]) / $n"
            for ii=max(1, i - nivel):min(ny, i + nivel)
                for jj=max(1, j - nivel):min(nx, j + nivel)
                    if !(ii == i && jj == j)
                        no_viz = nos[m[ii, jj]]
                        el_valido = true

                        el_tmp = Elemento(c, no_ref, no_viz, area, material)

                        for el in elementos
                            if tem_sobreposicao(el_tmp, el)
                                el_valido = false
                                break
                            end
                        end
                        if el_valido
                            push!(elementos, el_tmp)
                            c += 1
                        end
                    end
                end
            end
        end
    end
    return Estrutura(nos, elementos)
end

"""
    Otimizador(volume_max::Real, estrutura:: Estrutura)

Implementa as propriedades do modelo

volume_max: Volume de material disponível para o problema de otimização.
estrutura: Estrutura que será otimizada.
"""
mutable struct Otimizador
    volume_max::Float64
    estrutura::Estrutura
    x_min::Float64
    x_max::Float64
    x_ini::Float64
    move::Float64
    η::Float64
    x::Vector
    u::Vector

    function Otimizador(volume_max::Real, estrutura:: Estrutura)
        Ls = [comprimento(el) for el in estrutura.elementos]
        x_max = 1e3 * volume_max / minimum(Ls)
        x_ini = volume_max / sum(Ls)
        x_min = 0
        move = x_ini * 1e4
        η = 0.5
        x = fill(x_ini, length(Ls))
        u = Vector{Float64}(undef, num_gls_estrutura(estrutura))

        new(volume_max, estrutura, x_min, x_max, x_ini, move, η, x, u)
    end
end

"""
    fo(otimizador::Otimizador)::Float64

Retorna a função objetivo.
"""
function fo(otimizador::Otimizador)::Float64
    return vetor_forcas(otimizador.estrutura)' * otimizador.u
end

"""
    ∇fo(otimizador::Otimizador)::Float64

Retorna a sensibilidade da função objetivo.
"""
function ∇fo(otimizador::Otimizador)::Vector{Float64}
    els = otimizador.estrutura.elementos

    return [begin
        uₑ = otimizador.u[gls(els[i])]
        -uₑ' * k(els[i]) * uₑ
    end for i=eachindex(els)]
end

"""
    oc!(otimizador::Otimizador)

Atualiza as variáveis de projeto com o OC.
"""
function oc!(otimizador::Otimizador)
    n = length(otimizador.estrutura.elementos)

    sens_vol = [comprimento(el) for el in otimizador.estrutura.elementos]
    sens_fo = ∇fo(otimizador)
    vol_atual = volume(otimizador.estrutura)

    η = otimizador.η
    bm = -sens_fo ./ sens_vol

    l₁ = 0
    l₂ = 1.2 * maximum(bm)
    x_novo = Array{Float64}(undef, n)

    c = 1
    while (l₂ - l₁) > 1e-10 * (l₂ + 1) && c ≤ 1000
        lₘ = (l₂ + l₁) / 2

        # Correção de erro numérico
        be = max.(0, bm / lₘ)
        t = @. otimizador.x_min + (otimizador.x - otimizador.x_min) * be^η
        x_novo = @fastmath @. max(max(min(min(t, otimizador.x + otimizador.move), otimizador.x_max),
                                  otimizador.x - otimizador.move), otimizador.x_min)

        if (vol_atual - otimizador.volume_max) + sens_vol' * (x_novo - otimizador.x) > 0
            l₁ = lₘ
        else
            l₂ = lₘ
        end
        c += 1
    end
    otimizador.x = x_novo
end

"""
    otimizar!(otimizador::Otimizador; max_it=1000, ϵ=1e-9)

Minimiza a compliance do problema.
"""
function otimizar!(otimizador::Otimizador; max_it=1000, ϵ=1e-9)
    num_elems = length(otimizador.estrutura.elementos)
    erro = Inf

    iter = 0
    while erro > ϵ && iter ≤ max_it
        norm₀ = norm(otimizador.x)

        for i=1:num_elems
            otimizador.estrutura.elementos[i].area = otimizador.x[i]
        end

        otimizador.u = deslocamentos(otimizador.estrutura)

        oc!(otimizador)
        norm₁ = norm(otimizador.x)

        if iter > 1
            erro = abs(norm₁ - norm₀) / norm₀
        end

        iter += 1
        @info "It: $iter \tCompliance: $(fo(otimizador))\t Erro: $erro \tVolume: $(volume(otimizador.estrutura))"
    end
end

end

max_stars_repo_name: lorranfoliveira/ge

max_issues_repo_name: lorranfoliveira/ge

max_forks_repo_name: lorranfoliveira/ge


text:
using GeometryTypes, StaticArrays, ModernGL, Visualize
import GLAbstraction, GLWindow, ColorVectorSpace
import Transpiler: mix, smoothstep, gli
using Visualize: orthographicprojection, perspectiveprojection, lookat, normalmesh

function blinnphong2{NV, T}(V, N, L::Vec{NV, T}, color, ambient, light)
    diff_coeff = 2.3f0
    specular_power = 3.0f0
    specular = Vec3f0(0, 2, 3)
    spec_coeff = 0.77f0

    return ambient .* ambient .+
        color .* 1f0 .* color * diff_coeff .+
        specular .* specular_power .* specular * spec_coeff
end

# args = (Vec3f0(0), Vec3f0(0,1.5,3), Vec3f0(3), Vec3f0(1, 0, 0), shading, light)
# V, N, L, color, shading, light = args
# blinnphong2(args...)
# using Sugar, Transpiler
# m = Transpiler.GLMethod((blinnphong2, map(typeof, args)))
# ast = Sugar.getsource!(m)
# println(ast)
#
# ft = Sugar.expr_type(m, SlotNumber(13))
# args = map(x-> Sugar.expr_type(m, SSAValue(x)), 9:11)
# F = Sugar.instance(ft)
# typs = Tuple{ft, args...}
# m2 = Transpiler.GLMethod((broadcast, typs))
# ast2 = Sugar.sugared(m2.signature...)
# s2 = Sugar.getast!(m2)
# println(s2)
# Sugar.method_nargs(m2)
# @which Sugar.slotnames(m2)
# Sugar.getcodeinfo!(m2).slotnames
#
# F = Sugar.expr_type(m, SSAValue(0)).instance
#


text:
module MultiStochGrad

using Printf

using Logging
using Base.CoreLogging
using LinearAlgebra

debug_log = stdout
logger = ConsoleLogger(stdout, CoreLogging.Info)
global_logger(logger)

include("scsg.jl")
include("svrg.jl")

include("applis/LinearRegression.jl")
include("applis/LogisticRegression.jl")

# export list

export Observations,
    TermFunction,
    compute_value,
    TermGradient,
    compute_gradient!,
    Evaluator,
    get_nbterms,
    # from scsg
    SCSG,
    minimize,
    BatchSizeInfo,
    get_batchsizeinfo,
    #
    SVRG
end


text:
using BenchmarkTools
using StructuredProximalOperators
using Random

function main()

    benchGVal = true
    benchProx = true
    benchGrad = true
    benchHess = true

    g = regularizer_l1(2.5)
    n = 1000000
    Random.seed!(1567)
    x = rand(n).*5

    ## Prox
    y, M = prox_αg(g, x, 1.0)
    @show M

    if benchGVal
        println("\n + g(y)")
        g(y)
        @btime $g($y)
    end

    if benchProx
        println("\n + Prox computation")
        println("   - inplace")
        res = zeros(n)
        prox_αg!(g, res, x, 1.0)
        @btime prox_αg!($g, $res, $x, 1.0)

        println("   - with allocation")
        prox_αg(g, x, 1.0)
        @btime prox_αg($g, $x, 1.0)
    end

    ## Gradient, hessian
    if benchGrad
        println("\n + Gradient computation")
        println("   - inplace")
        res = zeros(n)
        ∇M_g!(g, M, res, x)
        @btime ∇M_g!($g, $M, $res, $x);

        println("   - with allocation")
        ∇M_g(g, M, x)
        @btime ∇M_g($g, $M, $x);
    end

    ## Gradient / hessian euclidean to riemannian
    if benchHess
        println("\n + Hessian-vector computation")
        println("   - inplace")
        d = zeros(n)
        res = zeros(n)
        ∇²M_g_ξ!(g, M, res, x, d)
        @btime ∇²M_g_ξ!($g, $M, $res, $x, $d);

        println("   - with allocation")
        ∇²M_g_ξ(g, M, x, d)
        @btime ∇²M_g_ξ!($g, $M, $res, $x, $d);
    end

    return
end

main()