text:
#Tile focused function prepare, and sort a tile etc.
"""
The grid_sale is defined as the twice the maximum of sum of theta and phi arm
The factor of 2 is used so that we can also obtain look at all potential collision
otherwise this factor of 2 is not needed for only to determine all possible assignments
"""
function get_grid_scale(fp_dic)
    grid_scale=2.0*maximum(fp_dic["LENGTH_R1"].+fp_dic["LENGTH_R2"])
    return grid_scale
end

"""sorts the targets by zone and then priority
"""
function sort_targets_ZonePriority!(targets_dic)
    #get the zone
    zone_xy=zeros(eltype(targets_dic["zone_x"]),targets_dic["ntarget"])
    get_zonexy!(targets_dic["zone_x"],targets_dic["zone_y"],targets_dic["fp_bound"],zone_xy)

    #Actual sorting, remapping and zone_map consistinting of indices range
    sort_ZoneXp_Apply!(zone_xy,targets_dic["PRIORITY"],targets_dic)

end

function sort_hardware_ZonePriority!(fp_dic,fp_bound)
    #get the zone
    nfp=size(fp_dic["zone_x"])
    zone_xy=zeros(eltype(fp_dic["zone_x"]),nfp)
    get_zonexy!(fp_dic["zone_x"],fp_dic["zone_y"],fp_bound,zone_xy)

    #Actual sorting, remapping and zone_map consistinting of indices range
    device_code=zeros(Int8,nfp)
    #POS is 0 and everything else is one as we only want to keep all positionr together
    #at the end and since it does reverse sort on xp
    device_code[findall(fp_dic["DEVICE_TYPE"] .!= "POS")] .= 1

    sort_ZoneXp_Apply!(zone_xy,device_code,fp_dic)
end

"""
sort the values based on combination of zone and XP
such that zones are together and XP is reverse ordered for each zone
This assumes zones are integers
It then apply the new sorting to everything in tdic which is array and of same size as Xp
Finally generate a map dictionary which when zone is specified gives the minimum and maximum index
One needs to be careful that zone_x and zone_y are not modified inside the function
as this function changes its arguments, one needs to make sure none of the Array is modified/sorted twice
even if they are refered within tdic and also provided directy
"""
function sort_ZoneXp_Apply!(zone_xy,Xp,tdic)

    #count number of entries
    nxp=size(Xp,1)

    #get the range of Xp
    xp_min=minimum(Xp)
    xp_max=maximum(Xp)

    #new quantity for sorting
    #The factor of 0.5 is to make sure the maximum value of priority is lower than next zone
    #It can be any value below 1.0, 0.5 is just a conveninet choice
    #This sort of decide the resolution of values in XP
    new_q=(Xp .- xp_min).*(0.5/(xp_max-xp_min))
    new_q=new_q+zone_xy
    izone_sort=sortperm(new_q,rev=:true)

    modified_ptr=[] #keep track of al the arrays sorted already to not apply sorting twice
    #Now apply this to all of the target arrays
    for tkey in keys(tdic)
        if(tdic[tkey] isa Union{BitArray,Array})
            #this_ptr=pointer_from_objref(tdic[tkey])
            this_ptr=repr(UInt64(pointer_from_objref(tdic[tkey])))
            if((size(tdic[tkey],1)==nxp) & (!(this_ptr in modified_ptr)) )
	        if(tdic[tkey] isa Union{BitMatrix,Matrix} )
	            tdic[tkey] .= tdic[tkey][izone_sort,:]
		else
	            tdic[tkey] .= tdic[tkey][izone_sort]
		end
		append!(modified_ptr,[this_ptr])
                #println("sorting $(tkey) $(this_ptr)")
	    elseif(!(tkey in ["fp_bound","file_map","ntarget"]))
                println("Warning: not sorting? ",tkey,size(tdic[tkey],1), this_ptr)
            end
        end
    end

    #Now generate the map from of indices for each zone
    #this_ptr=pointer_from_objref(zone_xy)
    this_ptr=repr(UInt64(pointer_from_objref(zone_xy)))
    if(!(this_ptr in modified_ptr) )
        zone_xy .= zone_xy[izone_sort]
        append!(modified_ptr,[this_ptr])
	#println("sorting zone_xy $(this_ptr)")
    end

    uszone=sort(unique(zone_xy),rev=false)
    first_zone=indexin(uszone,zone_xy)
    last_zone=(nxp+1) .- indexin(uszone,reverse(zone_xy))

    #generate the zone index map
    nunique_zone=size(uszone,1)

    tdic["zone_map"]=Dict()
    uu=1
    tdic["zone_map"][uszone[uu]]=[first_zone[uu],last_zone[uu]]
    max_per_zone=tdic["zone_map"][uszone[uu]][2]-tdic["zone_map"][uszone[uu]][1]
    min_per_zone=max_per_zone
    for uu in 2:nunique_zone
        tdic["zone_map"][uszone[uu]]=[first_zone[uu],last_zone[uu]]
        this_zone=last_zone[uu]-first_zone[uu]+1
        if(this_zone<min_per_zone)
            min_per_zone=this_zone
        elseif(this_zone>max_per_zone)
            max_per_zone=this_zone
        end
    end

    tdic["zone_map_keys"]=keys(tdic["zone_map"])
    tdic["count_maxper_zone"]=max_per_zone
    tdic["count_minper_zone"]=min_per_zone
end

function get_zonexy!(zone_x,zone_y,fp_bound,zone_xy)
    @. zone_xy=zone_x+ (zone_y*fp_bound.nx)
end

function get_zonexy(zone_x,zone_y,fp_bound)
    zone_xy=zone_x+ (zone_y*fp_bound.nx)
    return zone_xy
end

function split_zonexy(zone_xy,fp_bound)
    zone_x=zone_xy%fp_bound.nx
    zone_y=Int((zone_xy-zone_x)/fp_bound.nx)
    return zone_x,zone_y
end

"""
Loads all the tracer in a tile and attach a target_identification column to it
"""
function load_targets_intile(config,tile_id,tile_pass,group;
      columns=["X_UNIT","Y_UNIT","Z_UNIT","PRIORITY"],num_fba=1)
    targets_dic=Dict()
    targets_dic["file_map"]=Dict()

    col_accumulate=["index","filekey","collided"]
    for col in columns
        if(!(col in col_accumulate))
            append!(col_accumulate,[col])
        end
    end

    targets_dic["filekey"]=zeros(Int8,0)

    fno=1
    for tracer in keys(config["target"])
        #fname=config["target"][tracer]["JLDfile" ]
	#fba_jldfile=config["target"][tracer]["FBA_JLDfile" ]
	#fba_jldfile=TileFBA_FileName(config,tile_id)
        targets_dic["file_map"][fno]=[tracer,config["target"][tracer]["JLDfile" ]]

        tracer_tile=Load_tracers_intile(config,tracer,tile_id,tile_pass,group,columns;
			numobs_need=config["target"][tracer]["Num_obs"],num_fba=num_fba)

        nobj=size(tracer_tile["index"],1)
        for col in col_accumulate
            if(col=="filekey")
                targets_dic[col]=append!(targets_dic[col],(zeros(Int8,nobj) .+ fno))
            elseif(col in keys(targets_dic))
                #targets_dic[col]=append!(targets_dic[col],tracer_tile[col])
		targets_dic[col]=[targets_dic[col];tracer_tile[col]]
            else
                targets_dic[col]=tracer_tile[col]
            end
        end
        #println(tracer,' ',unique(targets_dic["PRIORITY"]),nobj,fno)
        fno +=1
	#println("load_targets_intile: ",tile_id,' ',num_fba,size(targets_dic["collided"]),
	#	' ',size(tracer_tile["collided"]))
    end
    targets_dic["ntarget"]=size(targets_dic["index"],1)
    return targets_dic
end

""" This prepare a tile for assignment given the config dictionary
This function execute the following steps
1) load targets in a given tile tile
"""
function Prepare_tile(config,tile_index,tiles_dic;group="",tile_date="2019-09-16T00:00:00",
        plate_scale=nothing,fp_dic=nothing,exc_dic=nothing,num_fba=1)
    #get the tile centre
    telra=tiles_dic["RA"][tile_index]
    teldec=tiles_dic["DEC"][tile_index]
    tile_id=tiles_dic["TILEID"][tile_index]
    tile_pass=tiles_dic["PASS"][tile_index]

    #load the plate_scale
    if(plate_scale==nothing)
        plate_scale=load_platescale(config["focal_plane"]["focalplane_dir"])
    end

    #Load the focal plane if needed
    if((fp_dic==nothing) || (exc_dic==nothing))
        fp_dic,exc_dic=load_hw_full_FocalPlane!(config["focal_plane"];date=tile_date)
    elseif((fp_dic["tile_date"]!= tile_date) || (exc_dic["tile_date"]!= tile_date) )
        fp_dic,exc_dic=load_hw_full_FocalPlane!(config["focal_plane"];date=tile_date)
    end

    #Load target
    targets_dic=load_targets_intile(config,tile_id,tile_pass,group;
                    columns=["X_UNIT","Y_UNIT","Z_UNIT","PRIORITY"],num_fba=num_fba)

    #convert to focal plane
    targets_dic["X_FP"], targets_dic["Y_FP"]= xyz_unit_2_xy_focalplane(telra, teldec,
    targets_dic["X_UNIT"],targets_dic["Y_UNIT"],targets_dic["Z_UNIT"],plate_scale)

    #define the grid scale
    grid_scale=get_grid_scale(fp_dic)

    #get the focal plane bound
    targets_dic["fp_bound"]=get_bound_focal_plane(targets_dic["X_FP"], targets_dic["Y_FP"],grid_scale;pad=1.0)

    #Assign zones
    targets_dic["zone_x"],targets_dic["zone_y"]=Assign_xy_zone(targets_dic["X_FP"], targets_dic["Y_FP"],targets_dic["fp_bound"])

    #Assifn zones for hardware in fp
    fp_dic["zone_x"],fp_dic["zone_y"]=Assign_xy_zone(fp_dic["OFFSET_X"], fp_dic["OFFSET_Y"],targets_dic["fp_bound"])

    #sort targets
    sort_targets_ZonePriority!(targets_dic)
    #sort devices
    sort_hardware_ZonePriority!(fp_dic,targets_dic["fp_bound"])

    return targets_dic,fp_dic,exc_dic
end

"""Runs the assignment for a given tile
num_fba: number of fba realization to generate
max_vertices maximum number of vertices any of the positioner polygon can have
max_ngb_pos: maximum number of negibohuring positioner any positioner can have
"""
function Assign_Tile(config,fp_dic,exc_dic,targets_dic;num_fba=3,max_vertices=30,max_ngb_pos=9)
    #Make sure the num_fba is propagated correctly
    @assert num_fba==size(targets_dic["collided"],2) "num_fba is not propagated correctl $(num_fba), $(size(targets_dic["collided"],2))"

    #generate two robot positioner to allocate memory to be used for efficiency
    pos_cur= POSRobotBody(max_vertices) #Assuming we will never have more than 30 vertices

    max_pos_zone=Integer(9)*fp_dic["count_maxper_zone"]
    pos_ngb_dic=Dict("index_central_positioner"=>0,"index"=>zeros(Int32,1max_pos_zone),
        "pos"=>[],"pos_index"=>[])
    for ii in 1:max_ngb_pos
        append!(pos_ngb_dic["pos"],[POSRobotBody(max_vertices)])
        append!(pos_ngb_dic["pos_index"],Int32(0))
    end

    #declare the strucutre to load for this particular zone with size of maximum possible
    #THis is so that only one time memory allocation is needed for efficiency
    #For any zone there can be 8 other neigbouring zone and hence a factor of 9
    max_target=Integer(9)*targets_dic["count_maxper_zone"]
    max_npriority=size(unique(targets_dic["PRIORITY"]),1)
    zone_target=Zone_Target(max_target,max_npriority)
    #to store the index values
    pot_ass_index=zeros(Int32,max_target)
    #pot_coll_index=zeros(Int32,max_target)
    #TO store either assigned or collided indx
    pot_ass_coll_index=zeros(Int32,max_target)
    # To store in_poly and in_circle behaviours
    in_poly=zeros(Bool,max_target)
    in_circle=zeros(Bool,max_target)

    #boolean for collision
    #collided=zeros(Bool,targets_dic["ntarget"],num_fba)
    #If something is already assigned needed number then it will be set to true
    #so that we do not assign this, for pass=0 this will be false everywhere in the begining
    collided = targets_dic["collided"]

    #integer for potential collision
    potential_coll=zeros(Bool,targets_dic["ntarget"],num_fba)

    #dictionary for collision locations
    #integer for Assignment
    position_assigned=zeros(Int32,size(fp_dic["OFFSET_X"],1),num_fba)
    if(config["OUTPUT"]["ThetaPhi_pos"])
        theta_assigned=zeros(Float64,size(fp_dic["OFFSET_X"],1),num_fba)
        phi_assigned=zeros(Float64,size(fp_dic["OFFSET_X"],1),num_fba)
    end

    local_ass=zeros(Int64,num_fba) #To keep the assignent in local index
    ind_ri=zeros(Int32,num_fba)

    #Pre-allocating some memory to avoid doing this millions of time for moving positioners
    move_diff_v=zeros(Float64,2)
    move_cos_sin=zeros(Float64,2)
    move_cos_sin_sum=zeros(Float64,2)
    move_cstheta=zeros(Float64,2)
    move_csphi=zeros(Float64,2)

    #for reproducibility we set the random seed
    Random.seed!(config["RandomSeed"]);

    for zone_xy in keys(fp_dic["zone_map"])
        #get the zone
        zone_x,zone_y=split_zonexy(zone_xy,targets_dic["fp_bound"])

        #load the targets in zones around
        load_target_zone(zone_target,zone_x,zone_y,targets_dic["fp_bound"],targets_dic)
        #unique sorted priority for this zone
        #us_priority=sort(unique(zone_target.PRIORITY[1:zone_target.ntarget[1]]),rev=true)

        #@show zone_target.ntarget,zone_target.usp

        tmp=0
        #loop over positioner in this zone
        for pp in fp_dic["zone_map"][zone_xy][1]:fp_dic["zone_map"][zone_xy][2]
            #load the info for this positioner
            #if((fp_dic["DEVICE_TYPE"][pp]!="POS") || (fp_dic["LENGTH_R1"][pp]==0.0) )
            if(!valid_positionter(fp_dic,pp;pos_type="assignment"))
                continue
            else #load positioner
                load_positioner!(pos_cur,pp,fp_dic,exc_dic)
                #load the negihbours of thie positioner
                load_positioner_neighbours!(pos_cur,pp,targets_dic["fp_bound"],fp_dic,exc_dic,pos_ngb_dic)
                #pos_x=fp_dic["OFFSET_X"][pp]
                #pos_y=fp_dic["OFFSET_Y"][pp]
            end

            #maximum distance to which positioner reach
            #max_r2=(fp_dic["LENGTH_R1"][pp]+fp_dic["LENGTH_R2"][pp])^2
            #minimum distance to which positioner can reach
            #min_r2=(fp_dic["LENGTH_R1"][pp]-fp_dic["LENGTH_R2"][pp])^2

            coll_r2=(sqrt(pos_cur.arms2_lim[2])+pos_cur.body_minmax[2])^2

            #scan for potential assignment and potential target
            for ii in 1:zone_target.ntarget[1]
                rad2=(zone_target.X_FP[ii]-pos_cur.xy[1])^2
                rad2 +=(zone_target.Y_FP[ii]-pos_cur.xy[2])^2
                if(rad2>coll_r2)#no collision
                    zone_target.pot_ass[ii]=false
                    zone_target.pot_coll[ii]=false
                elseif(rad2>pos_cur.arms2_lim[2]) #maximum arms length
                    zone_target.pot_ass[ii]=false
                    zone_target.pot_coll[ii]=true
                elseif(rad2>=pos_cur.arms2_lim[1]) #minimum arms length
                    zone_target.pot_ass[ii]=true
                    zone_target.pot_coll[ii]=false
                    #get theta_phi
                    zone_target.theta_pos[ii],zone_target.phi_pos[ii]=Hardware_xy_to_thetaphi(pos_cur,
                        zone_target.X_FP[ii],zone_target.Y_FP[ii],rad2)
                    #check if within theta phi range for the positioner
                    if(outside_theta_phi_range(zone_target.theta_pos[ii],zone_target.phi_pos[ii],
                                pos_cur.body_theta_limit,pos_cur.head_theta_limit) )
                        zone_target.pot_ass[ii]=false
                        zone_target.pot_coll[ii]=true
                    end
                else
                    zone_target.pot_ass[ii]=false
                    zone_target.pot_coll[ii]=false
                end
            end

            #a = @allocated begin
            #end; if a > 0 println("mem alloc: ",' ',a/1024/1024) end

            #priority_sel=Dict()
            #for up in us_priority
            #     priority_sel[up]=findall((zone_target.PRIORITY[1:zone_target.ntarget[1]].==up) .&
            #        zone_target.pot_ass[1:zone_target.ntarget[1]])
            #end
            #@show us_priority,priority_sel

            #possible_ass=zeros(Int64,100) #To store the actual possiblity

            #pot_ass_index=findall(zone_target.pot_ass[1:zone_target.ntarget[1]])
            #pot_coll_index=findall(zone_target.pot_coll[1:zone_target.ntarget[1]])
            npot_ass=find_index!(view(zone_target.pot_ass,1:zone_target.ntarget[1]),pot_ass_index)
            #npot_coll=find_index!(zone_target.pot_coll[1:zone_target.ntarget[1]],pot_coll_index)
            npot_ass_coll=find_index_OR!(view(zone_target.pot_ass,1:zone_target.ntarget[1]),
                view(zone_target.pot_coll,1:zone_target.ntarget[1]),pot_ass_coll_index)

            coll_map_dic=Dict()
            #loop over realizations
            for ri in 1:num_fba
                count=0
                for uu in 1:zone_target.ntarget[2] #do this in priority order high to low
                    #this_index=findall( zone_target.usp_bool_sel[1:zone_target.ntarget[1],uu] .&
                    #    zone_target.pot_ass[1:zone_target.ntarget[1]])
                    #println(' ',uu,' ',size(this_index),' ',
                    #    size(findall(zone_target.usp_bool_sel[1:zone_target.ntarget[1],uu])) )
                    for izone in pot_ass_index[1:npot_ass] #scan all the object in potential assignment
                        #@show pp,uu,ri,izone,zone_target.usp_bool_sel[izone,uu]
                        if(!zone_target.usp_bool_sel[izone,uu])
                            continue
                        end

                        #check if this object is collided
                        if(collided[zone_target.index[izone],ri])
                            continue
                        end

                        #@show ri,uu,izone
                        #check if this object is in potential collision
                        if(potential_coll[zone_target.index[izone],ri])
                            #Asuuming no-collision and should change if found collided
                            actual_collision=false

                            #Now check if this is actual collision with another positioner
                            #scan over neighbour positioner assigned
                            for ingb in 1:pos_ngb_dic["count_ngb"]
                                #check if positioner is assigned to a target
                                if(position_assigned[pos_ngb_dic["pos_index"][ingb],ri]==0)
                                    continue
                                end

                                #Target to which positioner is assigned
                                ingb_index=position_assigned[pos_ngb_dic["pos_index"][ingb],ri]
                                #check if there is an entry for this collision test
                                if( (izone,ingb,ingb_index) in keys(coll_map_dic) )
                                    actual_collision=coll_map_dic[(izone,ingb,ingb_index)]
                                else #now since it is not tested already check for this collision
                                    #move the neighbour to the assigned location
                                    if(config["OUTPUT"]["ThetaPhi_pos"])#use pre-computed values
                                        th_ngb=theta_assigned[pos_ngb_dic["pos_index"][ingb],ri]
                                        phi_ngb=phi_assigned[pos_ngb_dic["pos_index"][ingb],ri]
                                    else #compute the assigned theta_phi
                                        itr=position_assigned[pos_ngb_dic["pos_index"][ingb],ri]
                                        th_ngb,phi_ngb=Hardware_xy_to_thetaphi(pos_ngb_dic["pos"][ingb],
                                            targets_dic["X_FP"][itr],targets_dic["Y_FP"][itr])
                                    end
                                    #Now move the positioner
                                    move_positioner_to_theta_phi!(pos_ngb_dic["pos"][ingb],th_ngb,phi_ngb,
                                        diff_v=move_diff_v,cos_sin=move_cos_sin,cos_sin_sum=move_cos_sin_sum,
                                        cstheta=move_cstheta,csphi=move_csphi)

                                    #move current positioner to this target
                                    move_positioner_to_theta_phi!(pos_cur,
                                        zone_target.theta_pos[izone] ,zone_target.phi_pos[izone],
                                        diff_v=move_diff_v,cos_sin=move_cos_sin,cos_sin_sum=move_cos_sin_sum,
                                        cstheta=move_cstheta,csphi=move_csphi)

                                    #check for the collision
                                    actual_collision=positioner_collided(pos_cur,pos_ngb_dic["pos"][ingb])
                                    #store this in local map
                                    coll_map_dic[(izone,ingb,ingb_index)]=actual_collision

                                end

                                if(actual_collision)
                                    break
                                end
                            end #end of loop for scanning the neighbouring positioner

                            #check of the two positioners are collided
                            if(actual_collision)
                                continue
                            end
                        end

                        #This object can be assigned
                        count=count+1
                        #possible_ass[count]=izone
                        zone_target.possible_ass[count]=izone
                    end #end of izone loop

                    if(count>0)
                        #@show count,uu
                        izone_sel=zone_target.possible_ass[rand(1:count)]
                        local_ass[ri]=izone_sel
                        position_assigned[pp,ri]=zone_target.index[izone_sel]
                        if(config["OUTPUT"]["ThetaPhi_pos"])
                            theta_assigned[pp,ri]= zone_target.theta_pos[izone_sel]
                            phi_assigned[pp,ri]= zone_target.phi_pos[izone_sel]
                        end
                        break #to break the priority loop because we got assignment at this priority
                    else
                        local_ass[ri]=0
                        position_assigned[pp,ri]=0
                        #@show "nothing can be assigned",pp,ri,count
                    end
                end #end of uu loop
                #if(local_ass[ri]==0)
                #    @show "nothing can be assigned",pp,ri,count
                #end
            end #end of ri loop

            #break
            #udate the collision and potential collision
            uass=sort(unique(local_ass[local_ass .> 0]))

            #pot_ass_coll_index=sort([pot_ass_index[1:npot_ass];pot_coll_index[1:npot_coll]])
            #in_poly=zeros(Bool,size(pot_ass_coll_index,1))
            #in_circle=zeros(Bool,size(pot_ass_coll_index,1))
            for izone in uass
                #move positioner to this target
                move_positioner_to_theta_phi!(pos_cur,zone_target.theta_pos[izone],zone_target.phi_pos[izone],
                    diff_v=move_diff_v,cos_sin=move_cos_sin,cos_sin_sum=move_cos_sin_sum,
                    cstheta=move_cstheta,csphi=move_csphi)

                fill!(in_poly,false)
                fill!(in_circle,false)

                #determine the collisions and potential collision with body
                PointsInPoly_and_PaddedCircle!(view(zone_target.X_FP,pot_ass_coll_index[1:npot_ass_coll]),
                   view(zone_target.Y_FP,pot_ass_coll_index[1:npot_ass_coll]),
                    view(pos_cur.body_on_target,1:pos_cur.nseg[1],:),
                    pos_cur.head_minmax[2];max_rad=pos_cur.body_minmax[2],
                    in_poly=in_poly,in_circle=in_circle)

                #determine the collisions and potential collision with head
                PointsInPoly_and_PaddedCircle!(view(zone_target.X_FP,pot_ass_coll_index[1:npot_ass_coll]),
                   view(zone_target.Y_FP,pot_ass_coll_index[1:npot_ass_coll]),
                    view(pos_cur.head_on_target,1:pos_cur.nseg[2],:),
                    pos_cur.body_minmax[2];max_rad=pos_cur.head_minmax[2],
                    in_poly=in_poly,in_circle=in_circle)

                #ind_ri=findall(local_ass .== izone)
                nri=find_index!(local_ass .== izone,ind_ri)

                collided[zone_target.index[pot_ass_coll_index[in_poly]],ind_ri[1:nri]] .= true
                potential_coll[zone_target.index[pot_ass_coll_index[in_circle]],ind_ri[1:nri]] .= true
            end
            #@show position_assigned[pp,:],targets_dic["PRIORITY"][position_assigned[pp,:]]
            #tmp=tmp+1
            if(false)
                #plot_moved_positioner(pos_cur)
                #plot_zones(pos_cur,zone_target,fp_dic)
                ri=2
                nc=zone_target.ntarget[1]
                plot(zone_target.X_FP[1:nc],zone_target.Y_FP[1:nc],seriestype=:scatter,color=:black)

                plot!([targets_dic["X_FP"][position_assigned[pp,ri]]],
                    [targets_dic["Y_FP"][position_assigned[pp,ri]]],seriestype=:scatter,color=:red)

                plot_circle(pos_cur.xy[1],pos_cur.xy[2],sqrt(pos_cur.arms2_lim[2]),color=:red)
                plot_circle(pos_cur.xy[1],pos_cur.xy[2],sqrt(coll_r2),color=:blue)

                #plot the collided
                ipot_col=findall(potential_coll[:,ri])
                icol=findall(collided[:,ri])
                @show icol,targets_dic["X_FP"][icol]
                plot!(targets_dic["X_FP"][icol],targets_dic["Y_FP"][icol],seriestype=:scatter,color=:blue)
                plot!(targets_dic["X_FP"][ipot_col],targets_dic["Y_FP"][ipot_col],seriestype=:scatter,color=:green)
                #|> display

                move_positioner_to_theta_phi!(pos_cur,theta_assigned[pp,ri] ,phi_assigned[pp,ri],
                    diff_v=move_diff_v,cos_sin=move_cos_sin,cos_sin_sum=move_cos_sin_sum,
                    cstheta=move_cstheta,csphi=move_csphi)

                plot_moved_positioner(pos_cur;linecolor=:black,bodycolor=:blue,headcolor=:red,alpha=0.2,hold=true) |> display
            end

        end
        #break
    end

    if(config["OUTPUT"]["ThetaPhi_pos"])
        return position_assigned,theta_assigned,phi_assigned
    else
        return position_assigned,zeros(Float64,2),zeros(Float64,2)
    end
end

"""
takes the assignment, splits in different tracers and write a group in the JLD2 file for each tile
This is expected to be post-processed after each pass

first scans the file_map dictionary in targets_dic
for each key in the dictionary
    select all the object with non-zero assignments and write them to file
"""
function TileAssignment_To_JLD2File(tile_id,config,num_fba,fp_dic,targets_dic,position_assigned;
        sum_ass_global=zeros(Int16,2),ithis=nothing,target_assignment_global=nothing,
    location_assignment_global=nothing)

    ntarget=size(targets_dic["X_FP"],1)
    npos=size(position_assigned,1)

    #declare an array to transfer the assignment
    if(size(sum_ass_global,1)<ntarget)
        sum_ass=zeros(Int16,ntarget)
    else
        sum_ass=view(sum_ass_global,1:ntarget)
        fill!(sum_ass,0)
    end

    if(ithis==nothing)
        ithis=zeros(Int32,ntarget)
    end

    if(target_assignment_global==nothing)
        target_assignment=zeros(Bool,ntarget,num_fba)
        location_assignment=zeros(Int16,ntarget,num_fba) .- 1
    else
        target_assignment=view(target_assignment_global,1:ntarget,:)
        location_assignment=view(location_assignment_global,1:ntarget,:)
        fill!(target_assignment,false)
        fill!(location_assignment,-1)
    end

    for ri in 1:num_fba
        #println(ri,' ',maximum(position_assigned[:,ri]),' ',ntarget)
        @inbounds for ii in 1:npos
            if(position_assigned[ii,ri]>0)
                target_assignment[position_assigned[ii,ri],ri] = true
                location_assignment[position_assigned[ii,ri],ri]=fp_dic["LOCATION"][ii]
                sum_ass[position_assigned[ii,ri]] += 1
            end
        end
    end

    #histogram(sum_ass[sum_ass .>0],alpha=0.1)
    msg_out=""
    for tkey in keys(targets_dic["file_map"])
        tracer=targets_dic["file_map"][tkey][1]
        nthis=find_index_AND!((targets_dic["filekey"].==tkey),(sum_ass .>0),ithis)
        #println(tkey,targets_dic["file_map"][tkey],size(ithis,1))

        #histogram!(sum_ass[ithis],label=targets_dic["file_map"][tkey][1],alpha=0.3)

        #need to collect the assignment array and index array
        #index_arr=targets_dic["index"][ithis]
        #assignment_bool=target_assignment[ithis,:]

        #sort by index_arr
        isort=sortperm(targets_dic["index"][ithis[1:nthis]])
        #index_arr=index_arr[isort]
        #assignment_bool=assignment_bool[isort,:]

        #Now plot for testing
        #plot!(targets_dic["X_FP"][ithis][isort],targets_dic["Y_FP"][ithis][isort],
        #    seriestype=:scatter,label=targets_dic["file_map"][tkey][1],color=:auto,alpha=0.5) #|> display

        #get the filename to write
        fba_jldfile=TileFBA_FileName(config,tile_id)
	#fba_jldfile=config["target"][tracer]["FBA_JLDfile"]

        jldopen(fba_jldfile, "a+") do file
	   write(file,"$(tracer)_Bool_ass",target_assignment[ithis[isort],:])
	   write(file,"$(tracer)_location_ass",location_assignment[ithis[isort],:])
	   write(file,"$(tracer)_index",targets_dic["index"][ithis[isort]])
        end

	#number of assignment for the first realization
	msg_out=string(msg_out,tracer,":",sum(target_assignment[ithis[isort],1])," ")
        #println("written: ",fba_jldfile)
    end
    #println("sum: ",sum(sum_ass .==0),' ',minimum(sum_ass),' ',maximum(sum_ass))
    return msg_out
end

"""
Generate the file_name for a tile_output
"""
function TileFBA_FileName(config,tile_id)
    outdir=string(config["OUTPUT"]["JLD2_dir"],config["OUTPUT"]["FBA-tag"],"/")
    outfile=string(outdir,"TILEID_",tile_id,".jld2")
    return outfile
end

"""
file to write the fits output of each zone including the information about FBA
"""
function ZoneFITS_FileName(config,zone,tracer,group)
    outdir=string(config["OUTPUT"]["FITS_dir"],config["OUTPUT"]["FBA-tag"],"/tmp/")
    outfile=string(outdir,group,"_",tracer,"_zone_",zone,".fits.gz")
    return outfile
end

"""Runs the full assignmnet from disk-to-disk for a single tile
tile_index : index if the tile in tiles_dic to be assigned
verbose sets the verbosity of returned message
verbose=0 empty msg
verbos=1 overall time of the function
verbose>1 time for all steps
"""
function Run_Single_Tile(config,tile_index,tiles_dic;group="",tile_date="2019-09-16T00:00:00",
    plate_scale=nothing,fp_dic=nothing,exc_dic=nothing,verbose=0)

    events_dic=Dict()
    events_dic["events"]=[]
    events_dic["times"]=[]

    tile_id=tiles_dic["TILEID"][tile_index]
    if(verbose>0)
        append!(events_dic["times"],[now(UTC)])
        append!(events_dic["events"],["Time"])
    end

    #load data for the tile and focal plane
    targets_dic,fp_dic,exc_dic=Prepare_tile(config,tile_index,tiles_dic;group=group,tile_date=tile_date,
        plate_scale=plate_scale,fp_dic=fp_dic,exc_dic=exc_dic,num_fba=config["NumFBARealization"])
    if(verbose>1)
        append!(events_dic["times"],[now(UTC)])
        append!(events_dic["events"],["Prep"])
	#count the number of objects need assignment
	need_ass=targets_dic["ntarget"]-sum(targets_dic["collided"][:,1])
    end

    #println("Run Single: ",tile_id,' ',config["NumFBARealization"],size(targets_dic["collided"]))

    #run the assignment
    position_assigned,theta_assigned,phi_assigned=Assign_Tile(config,fp_dic,exc_dic,targets_dic,
        num_fba=config["NumFBARealization"],max_vertices=30,max_ngb_pos=9)

    if(verbose>1)
        append!(events_dic["times"],[now(UTC)])
        append!(events_dic["events"],["Ass"])
    end

    #write the data to file
    nass_str=TileAssignment_To_JLD2File(tile_id,config,config["NumFBARealization"],fp_dic,targets_dic,position_assigned)

    if(verbose>0)
        append!(events_dic["times"],[now(UTC)])
        if(verbose>1)
            append!(events_dic["events"],["Out"])
        end
    end

    msg_out=string(tile_id)
    #count the number of targets in tile vs number need assignment for the first realization
    if(verbose>1)
       msg_out=string(msg_out," targ:",need_ass,"/",targets_dic["ntarget"]," ",nass_str)
    end

    return convert_events_to_msg(events_dic,msg_out)
end

"""Runs the full assignmnet from disk-to-disk for a array of tiles
tile_index_arr : Array of indices if the tile in tiles_dic to be assigned
"""
function Run_Many_Tile(config,tile_index_arr,tiles_dic;group="",tile_date="2019-09-16T00:00:00",
    plate_scale=nothing,fp_dic=nothing,exc_dic=nothing,verbose=0,pre_msg="")

    #load plate_scale if needed
    if(plate_scale==nothing)
        plate_scale=load_platescale(config["focal_plane"]["focalplane_dir"])
    end

    #load the data for first tile to use things multiple times
    #load data for the tile and focal plane
    tile_index=tile_index_arr[1]
    targets_dic,fp_dic,exc_dic=Prepare_tile(config,tile_index,tiles_dic;group=group,tile_date=tile_date,
        plate_scale=plate_scale,fp_dic=fp_dic,exc_dic=exc_dic,num_fba=config["NumFBARealization"])

    #Now iterate over rest of the tile
    for tile_index in tile_index_arr[1:end]
        #output file name for this tile
        fba_jldfile=TileFBA_FileName(config,tiles_dic["TILEID"][tile_index])
	#check if this file exists then continue
	if(isfile(fba_jldfile))
	   println(pre_msg," TILEID:",tiles_dic["TILEID"][tile_index], " exists")
	   continue
	end

        #println(now(UTC)," Begin TILEID:",tiles_dic["TILEID"][tile_index])
        ret_msg=Run_Single_Tile(config,tile_index,tiles_dic;group=group,tile_date=tile_date,
        plate_scale=plate_scale,fp_dic=fp_dic,exc_dic=exc_dic,verbose=verbose)
        if(verbose>0)
            #println(now(UTC),pre_msg," TILEID:",ret_msg)
            println(pre_msg," TILEID:",ret_msg)
        end
    end
end

max_stars_repo_name: shadaba/MockFBA

max_issues_repo_name: shadaba/MockFBA

max_forks_repo_name: shadaba/MockFBA


text:
function get_sphere_uv(p::AbstractVec)
    ϕ::AbstractFloat = atan(z(p), x(p))
    Θ::AbstractFloat = asin(y(p))
    u::AbstractFloat = 1.0 - (ϕ + π) / 2π
    v::AbstractFloat = (Θ + 0.5π) / π
    u, v
end

struct Sphere <: Surface
    center::AbstractVec
    radius::AbstractFloat
end

function hit(sphere::Sphere, ray::Ray, range::NamedTuple{(:min, :max), Tuple{T, T}} where T <: AbstractFloat)::Tuple{Bool, AbstractFloat, Union{Nothing, HitPoint}}
    oc::AbstractVec = ray.origin - sphere.center
    a::AbstractFloat = dot(ray.direction, ray.direction)
    b::AbstractFloat = dot(oc, ray.direction)
    c::AbstractFloat = dot(oc, oc) - sphere.radius * sphere.radius
    discriminant::AbstractFloat = b * b - a * c

    if discriminant > 0
        t::AbstractFloat = (-b + copysign(sqrt(discriminant), -c)) / a;
        if range.min < t < range.max
            point::AbstractVec = point_at(ray, t)
            normal::AbstractVec = normalize(point - sphere.center)
            u::AbstractFloat, v::AbstractFloat = get_sphere_uv(normal)

            if c > 0
                return true, t, HitPoint(point, normal, u, v, true)
            else
                return true, t, HitPoint(point, -normal, u, v, false)
            end
        end
    end

    false, NaN, nothing
end

ext: jl

max_stars_repo_name: goooxu/JuliaRayTracing

max_issues_repo_name: goooxu/JuliaRayTracing

max_forks_repo_name: goooxu/JuliaRayTracing


text:
@testset "Widgets" begin

function Centered(layout, children...; props...)
    Box(layout, children...; valign = Center, halign = Center, props...)
end

function Item(name::String, children; centered::Bool = true, desc=name)
    Frame("", margin = 10, tooltip_text = desc,
        Box(:v, spacing = 10, margin = 10,
            (centered ? Centered(:v, margin=10, children) : children),
            Label(name, halign=Center),
        )
    )
end

function App(title; width, height)
    Window(title, width, height,
        Scrolled(
            Box(:v, spacing=5, margin = 20, valign=Center,
                Item("Label", desc = "Displays a small to medium amount of text",
                    Label("Lorem ipsum dolor...")
                ), # Item
                Item("ProgressBar", desc = "Displays progress visually" ,
                    ProgressBar(.50, hexpand=true)
                ), # Item
                Grid(spacing = 10,
                    Item("ImageView", desc = "Displays an image",
                        ImageView(asset"image.png")
                    ),
                    Item("Picture", desc = "Displays an image with custom width and height",
                        Picture(asset"image.png", 300, 300)
                    )
                ),
                Item("Slider", desc = "Horizontal slider",
                    centered=false,
                    Box(:h,
                        Slider(1:100, hexpand = true)
                    )
                ), # Item
                Grid(spacing = 10, halign = Center,
                    Item("Button", desc = "A button that emits a signal when clicked on",
                        Button("Button")
                    ), # Item
                    Item("ToggleButton", desc = "A button that maintain its state",
                        ToggleButton("ToggleButton")
                    ), # Item
                    Item("Dropdown", desc = "A widget that allows the user to choose an item from a list of options.",
                        Dropdown("Foo", "Bar", "Qux")
                    ), # Item
                    Item("ColorButton", desc = "Opens a color chooser dialog",
                        ColorButton("#f1a")
                    ), # Item
                    Item("Switch", desc = "A “light switch” style toggle",
                        Switch(false)
                    ) # Item
                ),
                Grid(spacing = 10, halign = Center,
                    Item("TextField", desc = "A single line text entry widget.",
                        TextField(text = "Hello, world!")
                    ), # Item
                    Item("SpinButton", desc = "Retrieves a number from the user",
                        SpinButton(-50:50)
                    ), # Item
                ),
                Frame("",
                    Expander("Expander", halign = Center, margin = 20,
                        Box(:v, spacing = 10, margin = 20,
                            Button("Foo"),
                            Button("Bar"),
                            Button("Qux"),
                        )
                    )
                ),
                Item("Grid",
                    Grid(spacing = 20,
                        Button("header") |> cspan(3),                         |,
                        Button("body")   |> span(2, 2), -, Button("items-a"), |,
                        -,                              -, Button("items-b"), |,
                        Button("footer") |> cspan(3),                         |,
                    )
                )
            ) # Box
        ) # Scrolled
    ) # Window
end

app = App("Widgets test", width=900, height=720)
@showall app
# @waitfor app.destroy
Lynx.destroy(app)

end # testset


text:
# subtract one unit to the current limit
# SI-
export psp_decrementilimit

const DECREMENT_ILIMIT = [0x53, 0x49, 0x2D, 0x0D]
"""
    psp_decrementilimit(io)

Subtract one unit to the current limit.
"""
function psp_decrementilimit(io_psp::IO)
  write(io_psp, DECREMENT_ILIMIT)
  return nothing
end


text:
# default values
force_scale = 200.0
tid = 400
fiber_size = 1
porder = 2
if length(ARGS) == 4
   global tid = parse(Int64, ARGS[1])
   global force_scale = parse(Float64, ARGS[2])
   global fiber_size = parse(Int64, ARGS[3])
   global porder = parse(Int64, ARGS[4])
end
printstyled("force_scale=$force_scale, tid=$tid, fiber_size=$fiber_size, porder=$porder\n", color=:green)
include("CommonFuncs.jl")
if Sys.MACHINE=="x86_64-apple-darwin18.6.0"
    matplotlib.use("macosx")
end
np = pyimport("numpy")

testtype = "PlaneStress"

#ρ = 1153.4 kg/m^3 = 1153.4 kg/m^3 * 0.000076 m = 0.0876584 kg/m^2
#E = 944836855.0kg/m/s^2 * 0.000076m = 71807.60098kg/s^2 = 0.07180760098kg/ms^2
prop = Dict("name"=> testtype, "rho"=> 0.0876584, "E"=>0.07180760098, "nu"=>0.4)

NT = Int64(1e5)
dt = 2.0e-3  #ms
T = NT * dt

nxc, nyc = 10,5
nx, ny =  nxc*fiber_size, nyc*fiber_size
Lx, Ly = 0.2, 0.1 #m
nodes, EBC, g, gt, FBC, fext, ft, npoints, node_to_point = BoundaryCondition(tid, nx, ny, porder, Lx, Ly; force_scale=force_scale)
elements = []
for j = 1:ny
    for i = 1:nx
        n = (nx*porder+1)*(j-1)*porder + (i-1)porder+1
        #element (i,j)
        if porder == 1
            #   4 ---- 3
            #
            #   1 ---- 2

            elnodes = [n, n + 1, n + 1 + (nx + 1), n + (nx + 1)]
        elseif porder == 2
            #   4 --7-- 3
            #   8   9   6
            #   1 --5-- 2
            elnodes = [n, n + 2, n + 2 + 2*(2*nx+1),  n + 2*(2*nx+1), n+1, n + 2 + (2*nx+1), n + 1 + 2*(2*nx+1), n + (2*nx+1), n+1+(2*nx+1)]
        else
            error("polynomial order error, porder= ", porder)
        end

        coords = nodes[elnodes,:]
        push!(elements,FiniteStrainContinuum(coords,elnodes, prop, 3))
    end
end

ndofs = 2
domain = Domain(nodes, elements, ndofs, EBC, g, FBC, fext)
setGeometryPoints!(domain, npoints, node_to_point)

state = zeros(domain.neqs)
∂u = zeros(domain.neqs)
globdat = GlobalData(state,zeros(domain.neqs), zeros(domain.neqs),∂u, domain.neqs, gt, ft)

assembleMassMatrix!(globdat, domain)
updateStates!(domain, globdat)

@show norm(domain.fext)

Δt = T/NT

####
# M = globdat.M
# _, stiff = assembleStiffAndForce(globdat, domain, 0.0)
# n = size(stiff)[1]
# M_d, stiff_d = Array(M), Array(stiff)
# # H1 q^{n+1} = H0 q^n + b
# # q = (dot{u}, u)

# H1 = zeros(Float64, 2*n, 2*n)
# H1[1:n,1:n] = M_d
# H1[1:n,n+1:2*n] = Δt/2.0*stiff_d
# H1[n+1:2*n,n+1:2*n] = M_d

# H0 = zeros(Float64, 2*n, 2*n)
# H0[1:n,1:n] = M_d
# H0[n+1:2*n,1:n] = Δt*M_d
# H0[1:n,n+1:2*n] = -Δt/2.0*stiff_d
# H0[n+1:2*n,n+1:2*n] = M_d - Δt^2/2.0*stiff_d

# b = zeros(Float64, 2*n)
# b[1:n] = Δt*domain.fext
# b[n+1:2*n] = Δt^2/2.0*domain.fext

# q = zeros(Float64, 2*n)
# q_err = zeros(Float64, 2*n)
##

#explicit solver
SolverInitial!(Δt, globdat, domain)

ω = EigenMode(Δt, globdat, domain)
@show "stable time step is ", 0.8 * 2/ω, " current time step is ", Δt
for i = 1:NT
    @info i, "/" , NT
    solver = ExplicitSolver(Δt, globdat, domain)
    if i%10 == 0
        ω = EigenMode(Δt, globdat, domain)
        @show "stable time step is ", 0.8 * 2/ω, " current time step is ", Δt
    end
end

# error()
if !isdir("$(@__DIR__)/Data")
    mkdir("$(@__DIR__)/Data")
end

if !isdir("$(@__DIR__)/Debug")
    mkdir("$(@__DIR__)/Debug")
end

if !isdir("$(@__DIR__)/Data/order$porder")
    mkdir("$(@__DIR__)/Data/order$porder")
end
if !isdir("$(@__DIR__)/Debug/order$porder")
    mkdir("$(@__DIR__)/Debug/order$porder")
end

write_data("$(@__DIR__)/Data/order$porder/$(tid)_$(force_scale)_$(fiber_size).dat", domain)
@save "Data/order$porder/domain$(tid)_$(force_scale)_$(fiber_size).jld2" domain

close("all")
frame = NT
# length scale(cm), stress scale(GPa), time scale(s)
scales = [100.0, 1.0/0.000076*1.0e-3, 1.0e-3]
visσ(domain, nx, ny,  domain.history["stress"][frame], domain.history["state"][frame]; scaling = scales)
xlabel("X (cm)")
ylabel("Y (cm)")
PyPlot.tight_layout()
savefig("Debug/order$porder/terminal$(tid)_$(force_scale)_$(fiber_size).png")

close("all")
Disp = zeros(Float64, NT+1, 5)
Disp[:,1] = LinRange(0, T, NT+1)
ux = [reshape(domain.history["state"][i][1:(nx*porder+1)*(ny*porder+1)], nx*porder+1, ny*porder+1)[end,end] for i = 1:length(domain.history["state"])]
uy = [reshape(domain.history["state"][i][(nx*porder+1)*(ny*porder+1)+1:end], nx*porder+1, ny*porder+1)[end,end] for i = 1:length(domain.history["state"])]
Disp[:,2], Disp[:,3] = ux, uy
ux = [reshape(domain.history["state"][i][1:(nx*porder+1)*(ny*porder+1)], nx*porder+1, ny*porder+1)[div((nx*porder+1),2),end] for i = 1:length(domain.history["state"])]
uy = [reshape(domain.history["state"][i][(nx*porder+1)*(ny*porder+1)+1:end], nx*porder+1, ny*porder+1)[div((nx*porder+1),2),end] for i = 1:length(domain.history["state"])]
Disp[:,4], Disp[:,5]  = ux, uy
plot(Disp[:, 1]*scales[3], Disp[:, 2]*scales[1], "--r",  label="Ref")
plot(Disp[:, 1]*scales[3], Disp[:, 3]*scales[1], "--y")
plot(Disp[:, 1]*scales[3], Disp[:, 4]*scales[1], "--b")
plot(Disp[:, 1]*scales[3], Disp[:, 5]*scales[1], "--g")

savefig("Debug/order$porder/disp$(tid)_$(force_scale)_$(fiber_size).png")


text:
"""
    GroupedFeatures(ps::AbstractVector{Int})

A type representing groups of features, wherein the first `ps[1]` features are one group,
    the next `ps[2]` features are the second group and so forth.
"""
struct GroupedFeatures{IV<:AbstractVector{Int}}
    ps::IV
    p::Int
    num_groups::Int
end

function GroupedFeatures(; group_size::Int, num_groups::Int)
    GroupedFeatures(fill(group_size, num_groups))
end

GroupedFeatures(ps) = GroupedFeatures(ps, sum(ps), length(ps))

ngroups(gr::GroupedFeatures) = gr.num_groups
nfeatures(gr::GroupedFeatures) = gr.p

function group_idx(gr::GroupedFeatures, i::Integer)
    starts = cumsum([1; gr.ps])[1:end-1]
    ends = cumsum(gr.ps)
    starts[i]:ends[i]
end

function group_summary(gr::GroupedFeatures, vec::AbstractVector, f)
    ps = gr.ps
    num_groups = gr.num_groups
    starts = cumsum([1; ps])[1:end-1]
    ends = cumsum(ps)
    el = eltype(f(vec))
    output = Vector{el}(undef, num_groups)
    for g = 1:num_groups
        output[g] = f(vec[starts[g]:ends[g]])
    end
    output
end

function group_expand(gr::GroupedFeatures, vec::AbstractVector)
    arr = zeros(eltype(vec), gr.p)
    for i = 1:gr.num_groups
        arr[group_idx(gr, i)] .= vec[i]
    end
    arr
end

function group_expand(gr::GroupedFeatures, el::Number)
    fill(el, gr.p)
end


text:
nodes = [1:29;]

adj = ones(29, 29) .* Inf
adj[1, [2, 3]] .= 1
adj[2, [1, 3, 4]] .= 1
adj[3, [1, 2, 4]] .= 1
adj[4, [2, 3, 5, 6, 7]] .= 1
adj[5, [4, 6]] .= 1
adj[6, [4, 5, 7, 9]] .= 1
adj[7, [4, 6, 8]] .= 1
adj[8, [7, 10]] .= 1
adj[9, [6, 10, 11]] .= 1
adj[10, [8, 9, 15]] .= 1
adj[11, [9, 12, 13, 15]] .= 1
adj[12, [11, 13, 18]] .= 1
adj[13, [11, 12, 14, 15, 18]] .= 1
adj[14, [13, 15, 16]] .= 1
adj[15, [10, 11, 13, 14, 16]] .= 1
adj[16, [14, 15, 17, 19, 21, 22]] .= 1
adj[17, [16, 18, 22]] .= 1
adj[18, [12, 13, 17, 23]] .= 1
adj[19, [16, 20, 21]] .= 1
adj[20, [19, 21, 26]] .= 1
adj[21, [16, 19, 20, 22, 25]] .= 1
adj[22, [16, 17, 21, 23, 25]] .= 1
adj[23, [18, 22, 24, 29]] .= 1
adj[24, [23, 25, 28]] .= 1
adj[25, [21, 22, 24, 26]] .= 1
adj[26, [20, 25, 27]] .= 1
adj[27, [26, 28]] .= 1
adj[28, [24, 27, 29]] .= 1
adj[29, [23, 28]] .= 1

types = Dict(
    1 => [8, 9, 13, 14, 24],
    2 => [
        1,
        2,
        3,
        4,
        5,
        6,
        7,
        10,
        11,
        12,
        15,
        16,
        17,
        18,
        19,
        20,
        21,
        22,
        23,
        25,
        26,
        27,
        28,
        29,
    ],
)


text:
# Skipping MacroDefinition: paMakeVersionNumber ( major , minor , subminor ) ( ( ( major ) & 0xFF ) << 16 | ( ( minor ) & 0xFF ) << 8 | ( ( subminor ) & 0xFF ) )
# Skipping MacroDefinition: paNoDevice ( ( PaDeviceIndex ) - 1 )
# Skipping MacroDefinition: paUseHostApiSpecificDeviceSpecification ( ( PaDeviceIndex ) - 2 )
# Skipping MacroDefinition: paFloat32 ( ( PaSampleFormat ) 0x00000001 )
# Skipping MacroDefinition: paInt32 ( ( PaSampleFormat ) 0x00000002 )
# Skipping MacroDefinition: paInt24 ( ( PaSampleFormat ) 0x00000004 )
# Skipping MacroDefinition: paInt16 ( ( PaSampleFormat ) 0x00000008 )
# Skipping MacroDefinition: paInt8 ( ( PaSampleFormat ) 0x00000010 )
# Skipping MacroDefinition: paUInt8 ( ( PaSampleFormat ) 0x00000020 )
# Skipping MacroDefinition: paCustomFormat ( ( PaSampleFormat ) 0x00010000 )
# Skipping MacroDefinition: paNonInterleaved ( ( PaSampleFormat ) 0x80000000 )
const paFloat32 = 0x00000001
const paNonInterleaved = 0x80000000

# Skipping MacroDefinition: paPlatformSpecificFlags ( ( PaStreamFlags ) 0xFFFF0000 )
# Skipping MacroDefinition: paInputUnderflow ( ( PaStreamCallbackFlags ) 0x00000001 )
# Skipping MacroDefinition: paInputOverflow ( ( PaStreamCallbackFlags ) 0x00000002 )
# Skipping MacroDefinition: paOutputUnderflow ( ( PaStreamCallbackFlags ) 0x00000004 )
# Skipping MacroDefinition: paOutputOverflow ( ( PaStreamCallbackFlags ) 0x00000008 )
# Skipping MacroDefinition: paPrimingOutput ( ( PaStreamCallbackFlags ) 0x00000010 )

const paNoFlag = 0x00000000
const paClipOff = 0x00000001
const paDitherOff = 0x00000002
const paNeverDropInput = 0x00000004
const paPrimeOutputBuffersUsingStreamCallback = 0x00000008

# Skipping Typedef: CXType_FunctionProto PaStreamCallback
const PaStreamCallback = Cvoid
# Skipping Typedef: CXType_FunctionProto PaStreamFinishedCallback
const PaStreamFinishedCallback = Cvoid


text:
# For comparison with ffr-LFDFT

push!(LOAD_PATH, "../LF_common/")
using m_LF3d
include("../LF_common/sparse_LF3d.jl")

const ⊗ = kron

function test_main( NN::Array{Int64} )

  println("NN = ", NN)

  Nx = NN[1]
  Ny = NN[2]
  Nz = NN[3]
  Npoints = Nx*Ny*Nz

  AA = [0.0, 0.0, 0.0]
  BB = [16.0, 16.0, 16.0]

  # Initialize LF
  LF = init_LF3d_c( NN, AA, BB, verbose=true )

  ∇2 = get_Laplacian3d_kron( LF )

  nzval = ∇2.nzval
  rowval = ∇2.rowval
  for ip = 1:size(nzval)[1]
    @printf("%5d %5d %18.10f\n", ip, rowval[ip], nzval[ip])
  end

  v = 1.1*ones(Npoints)
  Av = ∇2 * v

  for ip = 1:Npoints
    @printf("%5d %18.10f %18.10f\n", ip, v[ip], Av[ip])
  end

end

test_main( [3, 2, 4])


text:
function Ematrix(e, v, cond)
    if cond == "stress"
        E = (e / (1 - v * v)) * [
            1 v 0
            v 1 0
            0 0 (1 - v) / 2]
    elseif cond == "strain"
        E =
            (e * (1 - v) / ((1 + v) * (1 - 2 * v))) * [
                1 v / (v - 1) 0
                v / (v - 1) 1 0
                0 0 (1 - 2 * v) / (2 - 2 * v)]
    else
        println("no planestress or planestrain condition")
        exit()
    end
    return E
end
function quadra(uv,xs,ys)
    (u,v) = uv
    N1 = u*v*(u-1)*(v-1)/4; N5 = v*(1-u*u)*(v-1)/2;
    N2 = u*v*(u+1)*(v-1)/4; N6 = u*(u+1)*(1-v*v)/2;
    N3 = u*v*(u+1)*(v+1)/4; N7 = v*(1-u*u)*(v+1)/2;
    N4 = u*v*(u-1)*(v+1)/4; N8 = u*(u-1)*(1-v*v)/2;
            N9 = (1-u*u)*(1-v*v);
    x = dot(xs,[N1 N2 N3 N4 N5 N6 N7 N8 N9])
    y = dot(ys,[N1 N2 N3 N4 N5 N6 N7 N8 N9])
    return x,y
end
function Nuv(uv)
    (u,v) = uv
    Nu1 = (2*u-1)*(v-1)*v/4;  Nv1 = (u-1)*u*(2*v-1)/4;
    Nu2 = (2*u+1)*(v-1)*v/4;  Nv2 = (u+1)*u*(2*v-1)/4;
    Nu3 = (2*u+1)*(v+1)*v/4;  Nv3 = (u+1)*u*(2*v+1)/4;
    Nu4 = (2*u-1)*(v+1)*v/4;  Nv4 = (u-1)*u*(2*v+1)/4;
    Nu5 = -u*(v-1)*v;         Nv5 = -(u*u-1)*(2*v-1)/2;
    Nu6 = -(v*v-1)*(2*u+1)/2; Nv6 = -u*(u+1)*v;
    Nu7 = -u*(v+1)*v;         Nv7 = -(u*u-1)*(2*v+1)/2;
    Nu8 = -(v*v-1)*(2*u-1)/2; Nv8 = -u*(u-1)*v;
    Nu9 = 2*u*(v*v-1);        Nv9 = 2*v*(u*u-1);

    Nu = [Nu1 Nu2 Nu3 Nu4 Nu5 Nu6 Nu7 Nu8 Nu9]
    Nv = [Nv1 Nv2 Nv3 Nv4 Nv5 Nv6 Nv7 Nv8 Nv9]

    return [Nu; Nv]
end
function Jacobian(uv,xs,ys)
    J = Nuv(uv)*[xs' ys']
    return J
end
function Nxy(uv,xs,ys)
    J = Jacobian(uv,xs,ys)
    nuv = Nuv(uv)
    nxy = J\nuv
    return nxy
end
function N(uv)
    (u,v) = uv
    N1 = u*v*(u-1)*(v-1)/4; N5 = v*(1-u*u)*(v-1)/2;
    N2 = u*v*(u+1)*(v-1)/4; N6 = u*(u+1)*(1-v*v)/2;
    N3 = u*v*(u+1)*(v+1)/4; N7 = v*(1-u*u)*(v+1)/2;
    N4 = u*v*(u-1)*(v+1)/4; N8 = u*(u-1)*(1-v*v)/2;
            N9 = (1-u*u)*(1-v*v);

    N = [N1 N2 N3 N4 N5 N6 N7 N8 N9 zeros(1,9);
              zeros(1,9)           N1 N2 N3 N4 N5 N6 N7 N8 N9]
    return N
end
function LN(uv,xs,ys)
    nxy = Nxy(uv,xs,ys)
    nx = nxy[1,:]'
    ny = nxy[2,:]'
    nx = [nx zeros(1,9);
          zeros(1,9) nx]
    ny = [ny zeros(1,9);
          zeros(1,9) ny]

    LN = [1 0; 0 0; 0 1]*nx + [0 0; 0 1; 1 0]*ny

    return LN
end
function ke(E,xs,ys)
    g = [-0.949108 -0.741531 -0.405845 0.405845 0.741531 0.949108 0.0]
    w = [0.129485 0.279705 0.38183 0.38183 0.279705 0.129485 0.417959]
    kre = zeros(18,18)
    for i=1:7,j=1:7
        uv = [g[i] g[j]]
        ln = LN(uv,xs,ys)
        J = Jacobian(uv,xs,ys)
        kre = kre + w[i]*w[j]*ln'*E*ln*abs(det(J))
    end
    return kre
end
function IDs(bc)
    #note: dangerously RAMMY
    #returns
    #nbc: number of boundary conditions
    #ix: ids of edges for each condition
    #ni: contains number of edges for each condition
    nbc = count(x->x==-1,bc)+1
    ix = findall(x->x==-1,bc)
    ix = [ix[i][1] for i=1:nbc-1]
    prepend!(ix,0)
    ix = [collect(ix[i-1]+1:ix[i]-1) for i=2:nbc]

    ni = [length(i) for i in ix]
    ix = [ [i for i in bc[ix[j]]] for j=1:nbc-1]
    ix = [edge[ix[i],:] for i=1:nbc-1]
    return Int64(nbc-1),ni,ix
end
function BC(f,xedge,yedge,map)
    uv = [0. 0.]
    if map[1]==1 & map[2]==1
        ind = 1
        uv[2] = -1
    elseif map[2]==1 & map[3]==1
        ind = 2
        uv[1] = 1
    elseif map[3]==1 & map[4]==1
        ind = 1
        uv[2] = 1
    else
        ind = 2
        uv[1] = -1
    end
    g = [-0.949108 -0.741531 -0.405845 0.405845 0.741531 0.949108 0.0]
    w = [0.129485 0.279705 0.38183 0.38183 0.279705 0.129485 0.417959]
    rhs = zeros(18,1)
    for i=1:7
        uv[ind] = g[i]
        (x,y) = quadra(uv,xedge,yedge)
        nuv = Nuv(uv)
        xdot = nuv[ind,:]'*xedge'
        ydot = nuv[ind,:]'*yedge'
        road = (xdot[1]^2+ydot[1]^2)^(0.5)
        rhs = rhs + w[i]*N(uv)'*f([x y])*road
    end
    return rhs
end

max_stars_repo_name: billyon/first

max_issues_repo_name: billyon/first

max_forks_repo_name: billyon/first


text:
@testset "Compound node configuration tests" begin
    mpSim = MPsim()
    filePath = joinpath( "base", "baseConfig" )
    fileName = string( filePath, ".xlsx" )
    catFileName = joinpath( "base", "catalogue.xlsx" )

    XLSX.openxlsx( fileName ) do xf
        XLSX.openxlsx( catFileName ) do catXF
            MPIO.readAttributes( mpSim, xf, catXF )
            MPIO.readBaseNodes( mpSim, xf, catXF )
            MPIO.readCompoundNodes( mpSim, xf, catXF )
        end  # XLSX.openXLSX( catFileName ) do xf
    end  # XLSX.openxlsx( fileName ) do xf

    @test length( mpSim.compoundNodeList ) == 12
end  # @testset "Compound node configuration tests"


text:
# Solves sparse regression problems with the SLOPE method of Bogdan et al [Bogdan2014].
# Since the prox step in solving SLOPE is a quadratic program we here use Convex.jl
# to solve it. This should in general be an inferior solution to the FastProxSL1 method
# described in [Bogdan2014]. But it is fun and simple and still performs well. :)

# Beware! This is not yet bug-free...

# The key step in both the proximal gradient and the accelerated proximal gradient algorithms
# to solve SLOPE is to solve the prox problem. According to 2.2 and 2.3 in [Bogdan2014] the prox
# can be solved by a quadratic program which we implement in Convex.jl:
using Convex

# Set up x once and for all so we need not do it on each prox optimization round.
setup_x_and_constraints(p) = begin
  x = Variable(p)
  constraints = Constraint[]
  for i in 1:(p-1)
    constraints += (x[i] >= x[i+1])
  end
  constraints += (x[p] >= 0)
  return x, constraints
end

# The prox problem we solve assumes y are sorted and positive. We keep the permutation vector
# around so we can revert later.
sort_and_normalize(y) = begin
  absy = abs(y)
  return sortperm(absy, rev = true), sign(y), absy
end

# Go back from a sortperm to orig order.
function revert_sortperm(y, yperm)
  p = length(y)
  yr = zeros(p)
  for i in 1:p
    yi = yperm[i]
    yr[yi] = y[i]
  end
  yr
end

solve_prox_problem_with_convexjl(yproxvals, lambdas, x, constraints) = begin
  yperm, sngy, absy = sort_and_normalize(yproxvals)
  ynormalized = absy[yperm]

  proxproblem = minimize(0.5 * norm(ynormalized .- x, 2)^2 + sum(lambdas .* x), constraints)
  solve!(proxproblem)
  if proxproblem.status != :Optimal
    throw(string(proxproblem))
  end

  # Reorder and put back the signs...
  sngy .* revert_sortperm(x.value, yperm)
end

# The default choice for lambdas is Benjamini-Hochberg sequence where q is the FDR:
using Distributions
NormalDistr = Normal()
bh_lambdafunc(i, p, q = 0.10) = quantile(NormalDistr, 1 - i*q/(2*p))

# Page 7 of the latest Candes2015 paper on SLOPE states that a very simple lambda sequence can be used:
# But we found they need to be smaller for Convex.jl and FastProxSL1 prox solvers!?
candes2015_lambdafunc(k, p, q = 0.10) = sqrt(2*log(p/k))

# Use a given lambdas sequence or set it up based on the simple sequence from
# Candes2015 paper if not given.
function check_and_setup_lambdasequence(p, lambdas = nothing, q = 0.05)
  if lambdas != nothing
    @assert p == length(lambdas)
  else
    lambdas = map(k -> bh_lambdafunc(k, p, q), 1:p)
  end
  return lambdas
end

# The simplest SLOPE solver is just a proximal gradient. It assumes we have a way to solve the prox
# problem.
function solve_SLOPE_w_proximal_gradient(X::Matrix{Float64}, y::Vector{Float64}, proxsolver::Function;
  stepsizefunc = nothing,       # Default is to use 1/norm(X)^2, so a fixed scheme...
  lambdas = nothing,            # Default is to use the simple sequence below from Candes2015 paper...
  maxiterations = int(5e4), stopval = 1e-3, verbose = true)

  # Assert valid inputs, setup and precalc
  n, p = size(X)
  @assert n == length(y)
  lambdas = check_and_setup_lambdasequence(p, lambdas)
  Xprim = X'        # Precalc for speed
  bnew = bprev = randn(p)  # Initial guess is random
  if stepsizefunc == nothing
    limitval = (2 / norm(X)^2) / 2 # Divide by 2 to ensure it is less than the limitval
    stepsizefunc = (i,b) -> limitval
  end

  # Now iterate for maxiterations steps or until stopval reached.
  for k in 1:maxiterations
    tk = stepsizefunc(k, bprev)
    yproxvals = bprev .- tk * Xprim * (X * bprev .- y)
    bnew = proxsolver(yproxvals, lambdas)
    delta = norm(bnew .- bprev)
    println("k = ", k, ", delta = ", delta)
    if delta < stopval
      return bnew
    end
    bprev = bnew
  end
  return bnew
end

function solve_SLOPE_w_convexjl(X::Matrix{Float64}, y::Vector{Float64}; options...)
  p = size(X, 2)
  x, constraints = setup_x_and_constraints(p)
  proxsolver(yvals, lambdas) = solve_prox_problem_with_convexjl(yvals, lambdas, x, constraints)
  solve_SLOPE_w_proximal_gradient(X, y, proxsolver; options...)
end

# Now generate a random, sparse regression problem
function rand_sparse_problem(p; a = iceil(log(p)), n = int(p/2), amplitude = 10.0, sigma = amplitude/100, shuffle = true)
  beta = vcat(amplitude*ones(a)+randn(a), zeros(p-a))
  indices = collect(1:p)
  if shuffle
    shuffle!(indices)
    beta = beta[indices]
  end
  X = randn(n, p)
  errors = sigma*randn(n)
  y = X * beta .+ errors
  return beta, X, y, indices, n, p, a, sigma, amplitude, errors
end

beta, X, y, indices, n, p, a, sigma, amplitude, errors = rand_sparse_problem(100; shuffle = false);

# And let SLOPE loose on it...
#@time betahat = solve_SLOPE_w_convexjl(X, y)
#println(beta[1:a])
#println(betahat[1:a])
#println(norm(beta .- betahat))

# There are bugs though...

# Lets implement the stack-based FastProx SL1 alg of Bogdan et al to see if that helps...
# This is Alg 4 in the Bogdan2014 paper.
function fast_prox_sl1(yproxvals, lambdas)
  p = length(yproxvals)
  tuples = Any[]
  t = 0

  # Find optimal group levels
  for k in 1:p
    t += 1
    i, j, s = k, k, (yproxvals[k] - lambdas[k])
    w = abs(s)
    push!(tuples, (i, j, s, w))
    while (t > 1) && (tuples[t-1][4] .<= w)
      itm1, jtm1, stm1, wtm1 = tuples[t-1]
      it, jt, st, wt = tuples[t]
      wnew = (jtm1 - itm1 + 1)/(jt - itm1 + 1) * stm1 .+ (jt - it + 1)/(jt - itm1 + 1) * st
      tuples[t-1] = (itm1, jt, stm1 .+ st, abs(wnew))
      pop!(tuples) # Delete last entry
      t -= 1
    end
  end

  # Set entries in x for each block
  x = zeros(p)
  for l in 1:t
    for k in tuples[l][1]:tuples[l][2]
      x[k] = tuples[l][4]
    end
  end
  x
end

solve_prox_problem_with_fastproxsl1(yproxvals, lambdas) = begin
  yperm, sngy, absy = sort_and_normalize(yproxvals)
  ynormalized = absy[yperm]

  xvals = fast_prox_sl1(ynormalized, lambdas)

  # Reorder and put back the signs...
  sngy .* revert_sortperm(xvals, yperm)
end

function solve_SLOPE_w_fastproxsl1(X::Matrix{Float64}, y::Vector{Float64}; options...)
  solve_SLOPE_w_proximal_gradient(X, y, solve_prox_problem_with_fastproxsl1; options...)
end

# And let SLOPE loose on it...
@time betahat = solve_SLOPE_w_fastproxsl1(X, y)
println(betahat)
println(beta[1:a])
println(betahat[1:a])
println(length(find(abs(betahat) .> 0)))
println(norm(beta .- betahat))


text:
# This file is part of the TaylorIntegration.jl package; MIT licensed

# jetcoeffs!
doc"""
    jetcoeffs!(f, t, x)

Returns an updated `x` using the recursion relation of the
derivatives from the ODE $\dot{x}=dx/dt=f(t,x)$.

`f` is the function defining the RHS of the ODE, `x` is a `Taylor1{T}`
or a vector of that type, containing the Taylor expansion
of the dependent variables of the ODE and `t` is the independent
variable.
Initially, `x` contains only the 0-th order Taylor coefficients of
the current system state (the initial conditions), and `jetcoeffs!`
computes recursively the high-order derivates back into `x`.
"""
function jetcoeffs!{T<:Number}(eqsdiff, t0::T, x::Taylor1{T})
    order = x.order
    for ord in 1:order
        ordnext = ord+1

        # Set `xaux`, auxiliary Taylor1 variable to order `ord`
        @inbounds xaux = Taylor1( x.coeffs[1:ord] )

        # Equations of motion
        # TODO! define a macro to optimize the eqsdiff
        xdot = eqsdiff(t0, xaux)

        # Recursion relation
        @inbounds x.coeffs[ordnext] = xdot.coeffs[ord]/ord
    end
    nothing
end

function jetcoeffs!{T<:Number}(eqsdiff, t0::T, x::Vector{Taylor1{T}})
    order = x[1].order
    xaux = similar(x)
    for ord in 1:order
        ordnext = ord+1

        # Set `xaux`, auxiliary vector of Taylor1 to order `ord`
        @inbounds for j in eachindex(x)
            xaux[j] = Taylor1( x[j].coeffs[1:ord] )
        end

        # Equations of motion
        # TODO! define a macro to optimize the eqsdiff
        xdot = eqsdiff(t0, xaux)

        # Recursion relations
        @inbounds for j in eachindex(x)
            x[j].coeffs[ordnext] = xdot[j].coeffs[ord]/ord
        end
    end
    nothing
end

# stepsize
doc"""
    stepsize(x, epsilon)

Returns a time-step for a `x::Taylor1{T}` using a
prescribed absolute tolerance `epsilon`.
"""
function stepsize{T<:Number}(x::Taylor1{T}, epsilon::T)
    ord = x.order
    h = T(Inf)
    for k in (ord-1, ord)
        @inbounds aux = abs( x.coeffs[k+1] )
        aux == zero(T) && continue
        aux = epsilon / aux
        kinv = one(T)/k
        aux = aux^kinv
        h = min(h, aux)
    end
    return h
end

doc"""
    stepsize(q, epsilon)

Returns the minimum time-step for `q::Array{Taylor1{T},1}`,
using a prescribed absolute tolerance `epsilon`.
"""
function stepsize{T<:Number}(q::Array{Taylor1{T},1}, epsilon::T)
    h = T(Inf)
    for i in eachindex(q)
        @inbounds hi = stepsize( q[i], epsilon )
        h = min( h, hi )
    end
    return h
end

# taylorstep and taylorstep!
doc"""
    taylorstep(f, t0, x0, order, abs_tol)

Compute one-step Taylor integration for the ODE $\dot{x}=dx/dt=f(t, x)$
with initial condition $x(t_0)=x0$, returning the
time-step of the integration carried out and the updated value of `x0`.

Here, `x0` is the initial (and returned) dependent variable, `order`
is the degree used for the `Taylor1` polynomials during the integration
and `abs_tol` is the absolute tolerance used to determine the time step
of the integration.
"""
function taylorstep{T<:Number}(f, t0::T, x0::T, order::Int, abs_tol::T)
    # Initialize the Taylor1 expansions
    xT = Taylor1( x0, order )
    # Compute the Taylor coefficients
    jetcoeffs!(f, t0, xT)
    # Compute the step-size of the integration using `abs_tol`
    δt = stepsize(xT, abs_tol)
    x0 = evaluate(xT, δt)
    return δt, x0
end

doc"""
    taylorstep!(f, t0, x0, order, abs_tol)

Compute one-step Taylor integration for the ODE $\dot{x}=dx/dt=f(t, x)$
with initial conditions $x(t_0)=x0$, a vector of type T, returning the
step-size of the integration; `x0` is updated.

Here, `x0` is the initial (and updated) dependent variables, `order`
is the degree used for the `Taylor1` polynomials during the integration
and `abs_tol` is the absolute tolerance used to determine the time step
of the integration.
"""
function taylorstep!{T<:Number}(f, t0::T, x0::Array{T,1}, order::Int, abs_tol::T)
    # Initialize the vector of Taylor1 expansions
    xT = Array{Taylor1{T}}(length(x0))
    for i in eachindex(x0)
        @inbounds xT[i] = Taylor1( x0[i], order )
    end
    # Compute the Taylor coefficients
    jetcoeffs!(f, t0, xT)
    # Compute the step-size of the integration using `abs_tol`
    δt = stepsize(xT, abs_tol)
    evaluate!(xT, δt, x0)
    return δt
end

doc"""
    taylorstep(f, t0, t1, x0, order, abs_tol)

Compute one-step Taylor integration for the ODE $\dot{x}=dx/dt=f(t, x)$
with initial conditions $x(t_0)=x0$, returning the
time-step of the integration carried out and the updated value of `x0`.

Here, `x0` is the initial (and returned) dependent variables, `order`
is the degree used for the `Taylor1` polynomials during the integration
and `abs_tol` is the absolute tolerance used to determine the time step
of the integration. If the time step is larger than `t1-t0`, that difference
is used as the time step.
"""
function taylorstep{T<:Number}(f, t0::T, t1::T, x0::T, order::Int, abs_tol::T)
    @assert t1 > t0
    # Initialize the Taylor1 expansions
    xT = Taylor1( x0, order )
    # Compute the Taylor coefficients
    jetcoeffs!(f, t0, xT)
    # Compute the step-size of the integration using `abs_tol`
    δt = stepsize(xT, abs_tol)
    if δt ≥ t1-t0
        δt = t1-t0
    end
    x0 = evaluate(xT, δt)
    return δt, x0
end

doc"""
    taylorstep!(f, t0, t1, x0, order, abs_tol)

Compute one-step Taylor integration for the ODE $\dot{x}=dx/dt=f(t, x)$
with initial conditions $x(t_0)=x0$, a vector of type T, returning the
step-size of the integration carried out and updating `x0`.

Here, `x0` is the initial (and updated) dependent variables, `order`
is the degree used for the `Taylor1` polynomials during the integration
and `abs_tol` is the absolute tolerance used to determine the time step
of the integration. If the time step is larger than `t1-t0`, that difference
is used as the time step.
"""
function taylorstep!{T<:Number}(f, t0::T, t1::T, x0::Array{T,1},
        order::Int, abs_tol::T)
    @assert t1 > t0
    # Initialize the vector of Taylor1 expansions
    xT = Array{Taylor1{T}}(length(x0))
    for i in eachindex(x0)
        @inbounds xT[i] = Taylor1( x0[i], order )
    end
    # Compute the Taylor coefficients
    jetcoeffs!(f, t0, xT)
    # Compute the step-size of the integration using `abs_tol`
    δt = stepsize(xT, abs_tol)
    if δt ≥ t1-t0
        δt = t1-t0
    end
    evaluate!(xT, δt, x0)
    return δt
end

# taylorinteg
doc"""
    taylorinteg(f, x0, t0, tmax, order, abs_tol; keyword... )

This is a general-purpose Taylor integrator for the explicit ODE
$\dot{x}=f(t,x)$ with initial condition specified by `x0` at time `t0`.
It returns a vector with the values of time (independent variable),
and a vector (of type `typeof(x0)`) with the computed values of
the dependent variables. The integration stops when time
is larger than `tmax`, or the number of saved steps is larger
than `maxsteps`.

The integrator uses polynomial expansions on the independent variable
of order `order` and the parameter `abs_tol` serves to define the
time step using the last two Taylor coefficients of the expansions.

The current keyword argument is `maxsteps=500`.
"""
function taylorinteg{T<:Number}(f, x0::T, t0::T, t_max::T,
        order::Int, abs_tol::T; maxsteps::Int=500)

    # Allocation
    tv = Array{T}(maxsteps+1)
    xv = Array{T}(maxsteps+1)

    # Initial conditions
    nsteps = 1
    @inbounds tv[1] = t0
    @inbounds xv[1] = x0

    # Integration
    while t0 < t_max
        xold = x0
        δt, x0 = taylorstep(f, t0, x0, order, abs_tol)
        if t0+δt ≥ t_max
            x0 = xold
            δt, x0 = taylorstep(f, t0, t_max, x0, order, abs_tol)
            t0 = t_max
            nsteps += 1
            @inbounds tv[nsteps] = t0
            @inbounds xv[nsteps] = x0
            break
        end
        t0 += δt
        nsteps += 1
        @inbounds tv[nsteps] = t0
        @inbounds xv[nsteps] = x0
        if nsteps > maxsteps
            warn("""
            Maximum number of integration steps reached; exiting.
            """)
            break
        end
    end

    #return tv, xv
    return view(tv,1:nsteps), view(xv,1:nsteps)
end

function taylorinteg{T<:Number}(f, q0::Array{T,1}, t0::T, t_max::T,
        order::Int, abs_tol::T; maxsteps::Int=500)

    # Allocation
    tv = Array{T}(maxsteps+1)
    dof = length(q0)
    xv = Array{eltype(q0)}(dof, maxsteps+1)

    # Initial conditions
    @inbounds tv[1] = t0
    @inbounds xv[:,1] = q0[:]
    x0 = copy(q0)

    # Integration
    nsteps = 1
    while t0 < t_max
        xold = x0
        δt = taylorstep!(f, t0, x0, order, abs_tol)
        if t0+δt ≥ t_max
            x0 = xold
            δt = taylorstep!(f, t0, t_max, x0, order, abs_tol)
            t0 = t_max
            nsteps += 1
            @inbounds tv[nsteps] = t0
            @inbounds xv[:,nsteps] = x0[:]
            break
        end
        t0 += δt
        nsteps += 1
        @inbounds tv[nsteps] = t0
        @inbounds xv[:,nsteps] = x0[:]
        if nsteps > maxsteps
            warn("""
            Maximum number of integration steps reached; exiting.
            """)
            break
        end
    end

    #return tv, xv'
    return view(tv,1:nsteps), view(xv',1:nsteps,:)
end

# Integrate and return results evaluated at given time
function taylorinteg{T<:Number}(f, x0::T, trange::Range{T},
        order::Int, abs_tol::T; maxsteps::Int=500)

    # Allocation
    nn = length(trange)
    xv = Array{T,1}(nn)
    fill!(xv, T(NaN))

    # Initial conditions
    @inbounds xv[1] = x0

    # Integration
    iter = 1
    while iter < nn
        t0, t1 = trange[iter], trange[iter+1]
        nsteps = 0
        while nsteps < maxsteps
            xold = x0
            δt, x0 = taylorstep(f, t0, x0, order, abs_tol)
            if t0+δt ≥ t1
                x0 = xold
                δt, x0 = taylorstep(f, t0, t1, x0, order, abs_tol)
                t0 = t1
                break
            end
            t0 += δt
            nsteps += 1
        end
        if nsteps ≥ maxsteps && t0 != t1
            warn("""
            Maximum number of integration steps reached; exiting.
            """)
            break
        end
        iter += 1
        @inbounds xv[iter] = x0
    end

    return xv
end

function taylorinteg{T<:Number}(f, q0::Array{T,1}, trange::Range{T},
        order::Int, abs_tol::T; maxsteps::Int=500)

    # Allocation
    nn = length(trange)
    dof = length(q0)
    x0 = similar(q0, T, dof)
    fill!(x0, T(NaN))
    xv = Array{eltype(q0)}(dof, nn)
    @inbounds for ind in 1:nn
        xv[:,ind] = x0[:]
    end

    # Initial conditions
    @inbounds x0[:] = q0[:]
    @inbounds xv[:,1] = q0[:]

    # Integration
    iter = 1
    while iter < nn
        t0, t1 = trange[iter], trange[iter+1]
        nsteps = 0
        while nsteps < maxsteps
            xold = copy(x0)
            δt = taylorstep!(f, t0, x0, order, abs_tol)
            if t0+δt ≥ t1
                x0 = xold
                δt = taylorstep!(f, t0, t1, x0, order, abs_tol)
                t0 = t1
                break
            end
            t0 += δt
            nsteps += 1
        end
        if nsteps ≥ maxsteps && t0 != t1
            warn("""
            Maximum number of integration steps reached; exiting.
            """)
            break
        end
        iter += 1
        @inbounds xv[:,iter] = x0[:]
    end

    return xv'
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 538f05be-5053-11eb-19a0-959e34e1c2a1
using Pkg, DrWatson

# ╔═╡ 5d84f90c-5053-11eb-076b-5f30fc9685e3
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample, StanOptimize
	using StatisticalRethinking
	using StatsModelComparisons
end

# ╔═╡ 235c5298-5053-11eb-0608-435d1aa4716c
md" ## Clip-07-15s.jl"

# ╔═╡ 861accd4-5053-11eb-0432-81db212f4f38
begin
	sppnames = [:afarensis, :africanus, :hapilis, :boisei, :rudolfensis, :ergaster, :sapiens]
	brainvol = [438, 452, 612, 521, 752, 871, 1350]
	masskg = [37, 35.5, 34.5, 41.5, 55.5, 61, 53.5]
	df = DataFrame(species = sppnames, brain = brainvol, mass = masskg)
	scale!(df, [:mass, :brain])
	df.brain_std = df.brain/maximum(df.brain)
end;

# ╔═╡ 290d19ae-59cf-11eb-1636-772efccc7cb9
df

# ╔═╡ 24328f84-541d-11eb-0ec3-4df5b8a8cb19
begin
	scatter(df.mass, df.brain, xlab="body mass [kg]", ylab="brain vol [cc]",
		lab="Observations")
	for (ind, species) in pairs(df.species)
		annotate!([(df[ind, :mass] + 1, df[ind, :brain] + 30,
			Plots.text(df[ind, :species], 6, :red, :right))])
	end
	plot!()
end

# ╔═╡ ee275e7a-5067-11eb-325b-7760b758e85e
stan7_6 = "
data{
	int N;
	int K;
    vector[N] brain_std;
    matrix[N, K] mass;
}
parameters{
    real a;
    vector[K] b;
}
transformed parameters {
    vector[N] mu;
    mu = a + mass * b;
}
model{
    b ~ normal( 0 , 10 );
    a ~ normal( 0.5 , 1 );
    brain_std ~ normal( mu , 0.001 );
}
generated quantities{
    vector[N] log_lik;
    for (i in 1:N) log_lik[i] = normal_lpdf( brain_std[i] | mu[i], 0.001 );
}
";

# ╔═╡ 95cbb456-520e-11eb-2e0a-1d173ddc9dec
begin
	mass = create_observation_matrix(df.mass_s, 6)
	data = (N = 7, K = 6, brain_std = df.brain_std, mass = mass)
	m7_6s = SampleModel("m7.6s", stan7_6)
	rc7_6s = stan_sample(m7_6s; data=data)

	if success(rc7_6s)
		nt7_6s = read_samples(m7_6s, :namedtuple)
	end
end;

# ╔═╡ 172b171a-60b2-11eb-1d99-81f8e7cd6217
begin
	log_lik = nt7_6s.log_lik'
	n_sam, n_obs = size(log_lik)
	lppds = reshape(lppd(log_lik), n_obs)
end

# ╔═╡ 675109fc-59e6-11eb-247d-157746d5c626
sum(lppds)

# ╔═╡ e8343c86-5a79-11eb-1cd4-7b05ac85f7dd
md"
```
    mean   sd  5.5% 94.5% n_eff Rhat4
a   0.51 0.01  0.50  0.52    47  1.07
b1  0.88 0.01  0.86  0.90    50  1.07
b2  1.71 0.03  1.65  1.76    46  1.07
b3 -0.61 0.03 -0.65 -0.56    49  1.07
b4 -3.48 0.05 -3.56 -3.40    48  1.06
b5 -0.35 0.02 -0.38 -0.32    50  1.06
b6  1.63 0.02  1.59  1.66    49  1.06
```
"

# ╔═╡ bc6292b4-5a75-11eb-2a11-a31ae9ae5c4f
PRECIS(read_samples(m7_6s, :dataframe))

# ╔═╡ 77373916-5fe2-11eb-2df1-4bf002248937
begin
	loo, loos, pk = psisloo(log_lik)
	loo
end

# ╔═╡ c429b37a-5fe2-11eb-1345-190b929d58d8
sum(loos)

# ╔═╡ c429ee8a-5fe2-11eb-274e-31490ca32c5e
StatsModelComparisons.pk_qualify(pk)

# ╔═╡ c42aa528-5fe2-11eb-0492-1fc418032f9c
StatisticalRethinking.pk_plot(pk)

# ╔═╡ 7ce0cadc-61a2-11eb-1c20-353d2ad1df76
StatisticalRethinking.waic(log_lik)

# ╔═╡ bb780120-6ee5-11eb-2791-1903f6a16fa6
StatisticalRethinking.waic(m7_6s)

# ╔═╡ 91e9af0a-5065-11eb-212c-f751fd114263
md" ## End of clip-07-15s.jl"

# ╔═╡ Cell order:
# ╟─235c5298-5053-11eb-0608-435d1aa4716c
# ╠═538f05be-5053-11eb-19a0-959e34e1c2a1
# ╠═5d84f90c-5053-11eb-076b-5f30fc9685e3
# ╠═861accd4-5053-11eb-0432-81db212f4f38
# ╠═290d19ae-59cf-11eb-1636-772efccc7cb9
# ╠═24328f84-541d-11eb-0ec3-4df5b8a8cb19
# ╠═ee275e7a-5067-11eb-325b-7760b758e85e
# ╠═95cbb456-520e-11eb-2e0a-1d173ddc9dec
# ╠═172b171a-60b2-11eb-1d99-81f8e7cd6217
# ╠═675109fc-59e6-11eb-247d-157746d5c626
# ╟─e8343c86-5a79-11eb-1cd4-7b05ac85f7dd
# ╠═bc6292b4-5a75-11eb-2a11-a31ae9ae5c4f
# ╠═77373916-5fe2-11eb-2df1-4bf002248937
# ╠═c429b37a-5fe2-11eb-1345-190b929d58d8
# ╠═c429ee8a-5fe2-11eb-274e-31490ca32c5e
# ╠═c42aa528-5fe2-11eb-0492-1fc418032f9c
# ╠═7ce0cadc-61a2-11eb-1c20-353d2ad1df76
# ╠═bb780120-6ee5-11eb-2791-1903f6a16fa6
# ╟─91e9af0a-5065-11eb-212c-f751fd114263


text:
import Base64.base64encode
import Colors.hex

export to_svg

"""
    to_path_string(path)

Returns SVG path syntax for the given path.
"""
function to_path_string(path::Path)
    string("M", join(["$(point[1]),$(point[2])" for point in path], "L"))
end

"""
	to_svg(plot)

Returns an SVG representation of the given plot as a string.
"""
function to_svg(plot::PenPlot)
    ext = extent(plot)
    width = ext.lowerright[1] - ext.upperleft[1]
    height = ext.lowerright[2] - ext.upperleft[2]

    builder = IOBuffer()
    print(
        builder,
        """<svg
xmlns="http://www.w3.org/2000/svg"
height="400"
width="100%"
viewBox="$(ext.upperleft[1]) $(ext.upperleft[2]) $(width) $(height)">""",
    )

    for layer in plot.layers
        print(
            builder,
            """<g id="$(layer.name)"
stroke="#$(hex(layer.color))"
opacity="0.5"
fill="none">""",
        )

        for path in layer.paths
            print(
                builder,
                """<path
  d="$(to_path_string(path))"
  vector-effect="non-scaling-stroke"
  />""",
            )
        end

        print(builder, "</g>")
    end

    print(builder, "</svg>")

    String(take!(builder))
end

function Base.show(io::IO, mime::MIME"text/html", path::Path)
    plot = PenPlot([path])
    Base.show(io, mime, plot)
end

function Base.show(io::IO, mime::MIME"text/html", paths::MultiPath)
    plot = PenPlot(paths)
    Base.show(io, mime, plot)
end

function Base.show(io::IO, mime::MIME"text/html", plot::PenPlot)
    svgdata = to_svg(plot)

    # Write SVG element.
    write(io, svgdata)

    # Write download button.
    data_encoded = base64encode(svgdata)
    url = "data:image/svg+xml;base64,$(data_encoded)"
    write(io, "
")
    write(io, """<button>Download</button>""")
    write(io, "
")
end


text:
__precompile__()
module Unrolled

using MacroTools
using MacroTools: prewalk, postwalk

export @unroll, @code_unrolled
export unrolled_reduce, unrolled_filter, unrolled_intersect, unrolled_setdiff,
       unrolled_union, unrolled_in, unrolled_any, unrolled_all, unrolled_map,
       unrolled_foreach

function unrolled_filter end
""" `type_length(::Type)` returns the length of sequences of that type (only makes sense
for sequence-like types, obviously. """
function type_length end

"""
`type_size(::Type, dim)` returns the size of an array in the specified dimension.
"""
function type_size end

include("range.jl")

const expansion_funs = Dict{Function, Function}()

macro unroll_loop(niter_type::Type, loop)
    local niter
    @assert(@capture(loop, for var_ in seq_ loopbody__ end),
            "Internal error in @unroll_loop")
    try
        niter = type_length(niter_type)
    catch e
        # Don't unroll the loop, we can't figure out its length
        if isa(e, MethodError)
            return esc(loop)
        else rethrow() end
    end
    esc(:($Unrolled.@unroll_loop($niter, for $var in $seq
          $(loopbody...)
          end)))
end
macro unroll_loop(niter::Int, loop)
    @assert(@capture(loop, for var_ in seq_ loopbody__ end),
            "Internal error in @unroll_loop")
    esc(quote
        $([:(let $var = $seq[$i]; $(loopbody...) end) for i in 1:niter]...)
        nothing
        end)
end

macro unroll_loop(loop::Expr)
    @assert(@capture(loop, for var_ in 1:niter_ loopbody__ end),
            "Internal error in @unroll_loop")
    esc(quote $([:(let $var = $i; $(loopbody...) end) for i in 1:niter]...) end)
end

type_length(tup::Type{T}) where {T<:Tuple} = length(tup.parameters)
# Default fall-back
type_length(typ::Type) = length(typ)
type_size(typ::Type, i) = size(typ, i)

""" `function_argument_name(arg_expr)`

Returns the name (as a symbol) of this argument, where arg_expr is whatever can
be put in a function definition's argument list (eg. `len::Int=5`) """
function_argument_name(arg_expr) = MacroTools.splitarg(arg_expr)[1]

macro unroll(fundef)
    # This macro will turn the function definition into a generated function.
    di = splitdef(fundef)
    @assert !haskey(di, :params) "`@unroll` only supports parametric functions using `where ...` notation"
    fname = di[:name]
    args = di[:args]
    kwargs = get(di, :kwargs, [])
    body = di[:body]
    arg_vars = [a===nothing ? gensym() : a for a in map(function_argument_name, args)]
    kwarg_vars = map(function_argument_name, kwargs)
    all_args = [arg_vars; kwarg_vars]
    function seq_type(seq_var)
        @assert(seq_var in all_args,
                "Can only unroll a loop over one of the function's arguments")
        return Expr(:($), seq_var)
    end
    function seq_type_length(seq_var)
        @assert(seq_var in all_args,
                "Can only unroll a loop over one of the function's arguments")
        return Expr(:($), Expr(:call, :($Unrolled.type_length), seq_var))
    end
    function seq_type_size(seq_var, dim)
        @assert(seq_var in all_args,
                "Can only unroll a loop over one of the function's arguments")
        @assert(dim isa Integer,
                "Dimension argument must be an integer")
        return Expr(:($), Expr(:call, :($Unrolled.type_size), seq_var, dim))
    end
    process(x) = x
    function process(expr::Expr)
        if expr.args[1]==Symbol("@unroll")
            what = expr.args[3]
            @match what begin
                for var_ in 1:length(seq_) loopbody__ end =>
                    :($Unrolled.@unroll_loop(for $var in 1:$(seq_type_length(seq));
                                             $(loopbody...) end))
                for var_ in 1:size(seq_, dim_) loopbody__ end =>
                    :($Unrolled.@unroll_loop(for $var in 1:$(seq_type_size(seq, dim));
                                          $(loopbody...) end))
                for var_ in seq_ loopbody__ end =>
                    :($Unrolled.@unroll_loop($(seq_type(seq)),
                                             for $var in $seq; $(loopbody...) end))
                any_ => error("Cannot @unroll $what")
            end
        else
            expr
        end
    end
    # We walk over every expression in the function body, and replace the `@unroll`
    # loops with macros that will perform the actual unrolling (we use intermediate macros
    # for sanity)
    expansion = postwalk(process, body)
    exp_fun = Symbol(fname, :_unrolled_expansion_, gensym()) # gensym to support multiple methods
    return esc(quote
        # The expansion function (for easy calling)
        Base.@__doc__ function $exp_fun($(all_args...))
            $(Expr(:quote, expansion))
        end
        @generated function $fname($(args...); $(kwargs...)) where {$(di[:whereparams]...)}
            $exp_fun($(all_args...))
        end
        $Unrolled.expansion_funs[$fname] = $exp_fun
    end)
end

macro code_unrolled(expr)
    @assert(@capture(expr, f_(args__)))
    ar = gensym()
    esc(quote
        $ar = [$(args...)]
        macroexpand(@__MODULE__, $Unrolled.expansion_funs[$f](map(typeof, $ar)...))
    end)
end

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

@generated function unrolled_map(f, seq)
    :(tuple($((:(f(seq[$i])) for i in 1:type_length(seq))...)))
end
@generated function unrolled_map(f::F, seq1, seq2) where F
    @assert type_length(seq1) == type_length(seq2)
    :(tuple($((:(f(seq1[$i], seq2[$i])) for i in 1:type_length(seq1))...)))
end

@unroll function unrolled_foreach(f, seq)
    @unroll for x in seq; f(x) end
end

@generated function unrolled_reduce(f, v0, seq)
    niter = type_length(seq)
    expand(i) = i == 0 ? :v0 : :(f(seq[$i], $(expand(i-1))))
    return expand(niter)
end

function _unrolled_filter(f, tup)
    :($([Expr(:(...), :(f(tup[$i]) ? (tup[$i],) : ()))
         for i in 1:type_length(tup)]...),)
end
@generated unrolled_filter(f, tup) = _unrolled_filter(f, tup)
unrolled_intersect(tup1, tup2) = unrolled_filter(x->x in tup2, tup1)
unrolled_setdiff(tup1, tup2) = unrolled_filter(!(x->x in tup2), tup1)
unrolled_union() = ()
unrolled_union(tup1) = tup1
unrolled_union(tup1, tup2) = (tup1..., unrolled_setdiff(tup2, tup1)...)
unrolled_union(tup1, tup2, tupn...) =
    unrolled_reduce(unrolled_union, tup1, (tup2, tupn...))
""" `unrolled_in(obj, tup)` is like `in`. Beware that its return type is not
always known - see #21322 """
@inline @unroll function unrolled_in(obj, tup)
    @unroll for x in tup
        if obj == x
            return true
        end
    end
    return false
end

@unroll function unrolled_all(f, tup)
    @unroll for x in tup
        if !f(x)
            return false
        end
    end
    return true
end

@unroll function unrolled_any(f, tup)
    @unroll for x in tup
        if f(x)
            return true
        end
    end
    return false
end

end # module


text:
# Simplest test, just to check the most basic functionality.

using TrajOptSOCPs
using Test

@testset "Base Cases" begin
    # Write your tests here.
    @test fCheck2(4) == 16
    @test fCheck2(8) == 64
    @test fCheck2(10) == 100
end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
tol=0.00000001
b = [0.5]
c = [1.0, 0.3]

@assert minimum_value(b) - 0.5 < tol
@assert minimum_value(c) - 0.3 < tol

@assert algebraic_product(b) - 0.5 < tol
@assert algebraic_product(c) - 0.3 < tol

@assert bounded_difference(b) - 0.5 < tol
@assert bounded_difference(c) - 0.3 < tol

@assert drastic_product(b) - 0.5 < tol
@assert drastic_product(c) - 0.3 < tol

@assert einstein_product(b) - 0.5 < tol
@assert einstein_product(c) - 0.3 < tol

@assert hamacher_product(b) - 0.5 < tol
@assert hamacher_product(c) - 0.3 < tol


text:
function paint_X!(r::Robot)
    paint_X_Sud_right(r)
    paint_X_Sud_left(r)
    paint_X_Nord_left(r)
    paint_X_Nord_right(r)
    putmarker!(r)
end

function paint_X_Sud_right(r::Robot)
    while isborder(r, Nord) == false
        if isborder(r, Sud) == false
            if isborder(r, Ost) == false
                move!(r, Ost)
                move!(r, Sud)
                putmarker!(r)
            else
                break
            end
        else
            break
        end
    end

    for i in 1:100
        move!(r, Nord)
        move!(r, West)
        if ismarker(r) == false
            break
        end
    end
end

function paint_X_Sud_left(r::Robot)
    while isborder(r, Nord) == false
        if isborder(r, Sud) == false
            if isborder(r, West) == false
                move!(r, West)
                move!(r, Sud)
                putmarker!(r)
            else
                break
            end
        else
            break
        end
    end

    for i in 1:100
        move!(r, Nord)
        move!(r, Ost)
        if ismarker(r) == false
            break
        end
    end
end

function paint_X_Nord_left(r::Robot)
    while isborder(r, Sud) == false
        if isborder(r, Nord) == false
            if isborder(r, West) == false
                move!(r, West)
                move!(r, Nord)
                putmarker!(r)
            else
                break
            end
        else
            break
        end
    end

    for i in 1:100
        move!(r, Sud)
        move!(r, Ost)
        if ismarker(r) == false
            break
        end
    end
end

function paint_X_Nord_right(r::Robot)
    while isborder(r, Sud) == false
        if isborder(r, Nord) == false
            if isborder(r, Ost) == false
                move!(r, Ost)
                move!(r, Nord)
                putmarker!(r)
            else
                break
            end
        else
            break
        end
    end

    for i in 1:100
        move!(r, Sud)
        move!(r, West)
        if ismarker(r) == false
            break
        end
    end
end


text:
export FakeFracFldMat

################################################################################
#
#  FakeFracFldMat
#
################################################################################

mutable struct FakeFracFldMat
  num::MatElem
  den::RingElem
  rows::Int
  cols::Int

  function FakeFracFldMat()
    z = new()
    return z
  end

  function FakeFracFldMat(x::MatElem, y::RingElem, simplified::Bool = false)
    z = new()
    z.num = x
    z.den = y
    z.rows = nrows(x)
    z.cols = ncols(x)
    if !simplified
      simplify_content!(z)
    end
    return z
  end

  function FakeFracFldMat(x::Tuple{MatElem, RingElem}, simplified::Bool = false)
    z = new()
    z.num = x[1]
    z.den = x[2]
    z.rows = nrows(x[1])
    z.cols = ncols(x[1])
    if !simplified
      simplify_content!(z)
    end
    return z
  end

  # TODO: Maybe this should be a copy option
  function FakeFracFldMat(x::MatElem)
    z = new()
    z.num = x
    z.den = one(base_ring(x))
    z.rows = nrows(x)
    z.cols = ncols(x)
    return z
  end
end

################################################################################
#
#  FakeFracFldMat.jl : Model matrices of the fraction field of an integral domain R as matrices in R  with denominator in R
#
################################################################################

Base.iszero(x::FakeFracFldMat) = iszero(x.num)

function Base.numerator(x::FakeFracFldMat; copy::Bool = true)
  if copy
    return deepcopy(x.num)
  else
    return x.num
  end
end

function Base.denominator(x::FakeFracFldMat; copy::Bool = true)
  if copy
    return deepcopy(x.den)
  else
    return x.den
  end
end

Hecke.nrows(x::FakeFracFldMat) = nrows(x.num)

Hecke.ncols(x::FakeFracFldMat) = ncols(x.num)

function simplify_content!(x::FakeFracFldMat)
  c = content(x.num)
  c = gcd(c, x.den)
  if !isone(c)
    x.num = divexact(x.num, c)
    x.den = divexact(x.den, c)
  end
  y = canonical_unit(x.den)
  if !isone(y)
    mul!(x.den, x.den, y)
    mul!(x.num, x.num, y)
  end
end

################################################################################
#
#  Getindex
#
################################################################################

Base.getindex(a::FakeFracFldMat, i::Int, j::Int) = FractionField(base_ring(a))(a.num[i, j], a.den)

################################################################################
#
#  Unary operations
#
################################################################################

function Base.:(-)(x::FakeFracFldMat)
  return FakeFracFldMat(-x.num, x.den, true)
end

#TODO: may be possible to simplify more efficiently.
#The content of the numerator of the inverse may be non trivial!
function Base.inv(x::FakeFracFldMat)
  i, d_i = pseudo_inv(x.num)
  g = gcd(d_i, x.den)
  if isone(g)
    return FakeFracFldMat(i * x.den, d_i)
  elseif g == d_i
    return FakeFracFldMat(i * divexact(x.den, d_i))
  elseif g == x.den
    return FakeFracFldMat(i, divexact(d_i, x.den))
  else
    return FakeFracFldMat(i * divexact(x.den, g), divexact(d_i, g))
  end
end

################################################################################
#
#  Binary operations
#
################################################################################

function Base.:(+)(x::FakeFracFldMat, y::FakeFracFldMat)
  t = y.den * x.num + x.den * y.num
  d = x.den * y.den
  z = FakeFracFldMat(t, d)
  return z
end

function Base.:(*)(x::FakeFracFldMat, y::FakeFracFldMat)
  t = x.num * y.num
  d = x.den * y.den
  z = FakeFracFldMat(t, d)
  return z
end

function Hecke.mul!(z::FakeFracFldMat, x::FakeFracFldMat, y::FakeFracFldMat)
  @assert ncols(x) == nrows(y) "Wrong dimensions"
  @assert nrows(z) == nrows(x) "Wrong dimensions"
  @assert ncols(z) == ncols(y) "Wrong dimensions"
  mul!(z.num, x.num, y.num)
  mul!(z.den, x.den, y.den)
  simplify_content!(z)
  return z
end

################################################################################
#
#  Adhoc binary operations
#
################################################################################

function Base.:(*)(x::FakeFracFldMat, y::MatElem)
  t = x.num * y
  z = FakeFracFldMat(t, denominator(x))
  return z
end

function Base.:(*)(x::MatElem, y::FakeFracFldMat)
  t = x * y.num
  z = FakeFracFldMat(t, denominator(y))
  return z
end

function Base.:(*)(x::RingElem, y::FakeFracFldMat)
  g = gcd(x, y.den)
  if isone(g)
    return FakeFracFldMat(y.num * x, y.den, true)
  else
    return FakeFracFldMat(y.num * divexact(x, g), divexact(y.den, g), true)
  end
end

Base.:(*)(x::FakeFracFldMat, y::RingElem) = y * x

################################################################################
#
#  Comparison
#
################################################################################

Base.:(==)(x::FakeFracFldMat, y::FakeFracFldMat) = (x.num == y.num) && (x.den == y.den)

Base.isequal(x::FakeFracFldMat, y::FakeFracFldMat) = (x.num == y.num) && (x.den == y.den)

################################################################################
#
#  Determinant
#
################################################################################

function Hecke.det(x::FakeFracFldMat)
  nrows(x) != ncols(x) && error("Matrix must be square")

  return det(x.num)//(x.den)^nrows(x)
end

################################################################################
#
#  Rank
#
################################################################################

Hecke.rank(x::FakeFracFldMat) = rank(x.num)

################################################################################
#
#  Transpose
#
################################################################################

function Hecke.transpose(M::FakeFracFldMat)
  return FakeFracFldMat(transpose(numerator(M, copy = false)), denominator(M))
end


text:
include("constants.jl")

function init_config!()
    from_location = joinpath(PKG_ROOT, NEW_APP_PATH)
    dir_cp(from_location, PROJ_ROOT)
end

function dir_cp(src::String, dest::String)
    items = readdir(src)
    dirs = filter(item -> isdir(joinpath(src,item)), items)
    files = filter(item -> isfile(joinpath(src,item)), items)
    for dir in dirs
        newpath = joinpath(dest, dir)
        if !isdir(newpath)
            mkdir(newpath)
        end
    end
    for file in files
        destination = joinpath(dest, file)
        if !isfile(destination)
            cp(joinpath(src,file), destination)
        end
    end
    for dir in dirs
        dir_cp(joinpath(src, dir), joinpath(dest, dir))
    end
end


text:
########################################################################
#################### Summary Tables ####################################
########################################################################

count(x) = mapreduce(y->1, +, x)
p10(x) = quantile(x |> collect, .1)
p25(x) = quantile(x |> collect, .25)
p50(x) = quantile(x |> collect, .50)
p75(x) = quantile(x |> collect, .75)
p90(x) = quantile(x |> collect, .9)

function default_stats(detail::Bool)
    if !detail
        return ("Obs"         => count ∘ skipmissing,
                "Mean"        => mean ∘ skipmissing,
                "Std. Dev."   => std ∘ skipmissing,
                "Min"         => minimum ∘ skipmissing,
                "Max"         => maximum ∘ skipmissing)
    else
        return ("Obs"         => count ∘ skipmissing,
                "Mean"        => mean ∘ skipmissing,
                "Std. Dev."   => std ∘ skipmissing,
                "Min"         => minimum ∘ skipmissing,
                "p10"         => p10 ∘ skipmissing,
                "p25"         => p25 ∘ skipmissing,
                "p50"         => p50 ∘ skipmissing,
                "p75"         => p75 ∘ skipmissing,
                "p90"         => p90 ∘ skipmissing,
                "Max"         => maximum ∘ skipmissing)
    end
end

NumericCol = AbstractVector{T} where {T1<:Real, T2<:Real,
                                      T<:Union{T1, Union{T2, Missing}}}

tuplefy(x) = tuple(x)
tuplefy(x::Tuple) = x

function promotearray(x::AbstractArray{S, N}) where {S,N}
    types = unique(typeof(val) for val in x)
    T     = reduce(promote_type, types)
    return Array{T,N}
end

function summarize(df::AbstractDataFrame, fields=names(df);
                   detail=false, stats=default_stats(detail), kwargs...)

    # Determine whether each column is numeric or not
    numeric = Dict(header => typeof(df[!,header]) <: NumericCol
                   for header in fields)
    cols = TableCol[]
    for pair in tuplefy(stats)
        col = TableCol(pair.first)
        for header in fields
            if numeric[header]
                col[header] = pair.second(df[!,header])
            else
                col[header] = ""
            end
        end
        push!(cols, col)
    end

    return hcat(cols...)
end

function summarize(df::AbstractDataFrame, field::Symbol; kwargs...)
    summarize(df, vcat(field); kwargs...)
end

function summarize_by(df, byname::Symbol,
                      fields=setdiff(names(df), vcat(string.(byname)));
                      kwargs...)
    tabs = []
    gd = groupby(df, byname)
    for sub in gd
        tab = summarize(sub, fields; kwargs...)
        vals= unique(sub[!,byname])
        length(vals) == 1 || throw(error("Groupby isn't working"))
        idx = vals[1]
        push!(tabs, string(idx)=>tab)
    end

    return append_table(tabs...)
end

########################################################################
#################### Cross Tabulations #################################
########################################################################

function tabulate(df::AbstractDataFrame, field::Symbol)

    # Count the number of observations by `field`
    tab = combine(groupby(df, field), field => length => :_N)

    # Construct a Frequency Column
    sort!(tab, field)
    vals  = tab[!,field] .|> Symbol
    freq  = tab[!,:_N]
    pct   = freq/sum(freq)*100
    cum   = cumsum(pct)

    # Construct Table
    col1 = append_table(TableCol("Freq.",    vals, freq),
                        TableCol("Freq.", "Total"=>sum(freq)))
    col2 = append_table(TableCol("Percent",  vals, pct ),
                        TableCol("Percent", "Total"=>sum(pct)))
    col3 = TableCol("Cum.",     vals, cum )
    col = hcat(col1, col2, col3)
end

function tabulate(df::AbstractDataFrame, field1::Symbol, field2::Symbol)

    # Count the number of observations by `field`
    fields = vcat(field1, field2)
    df     = dropmissing(df[!,fields], disallowmissing=true)
    tab = combine(groupby(df, fields), field1 => length => :_N)
    sort!(tab, [field1, field2])

    # Put it into wide form
    tab = unstack(tab, field1, field2, :_N)

    # Construct the table
    vals = Symbol.(sort(unique(df[!,field2])))
    cols = []
    for val in vals
        col  = TableCol(val, Vector(tab[!,field1]), tab[!,val])
        col2 = TableCol(val, "Total" => sum(coalesce.(tab[!,val], 0)))
        push!(cols, append_table(field1=>col, ""=>col2))
    end

    sums = sum(coalesce.(Matrix(tab[!,vals]), 0), dims=2) |> vec
    tot1 = TableCol("Total", Vector(tab[!,field1]), sums)
    tot2 = TableCol("Total", "Total" => sum(sums))
    tot  = append_table(field1=>tot1, ""=>tot2)

    ret  = join_table(field2=>hcat(cols...), tot)
end