text:
function assemble_structures(layup,n_pt,twist,sloc,xaf,yaf,nchord,lam_t,designparams,x_offset,orientation,xloc,zloc,geom,webloc)
    # assemble structural properties
    mat = Array{Array{Composites.material,1}}(n_pt)
    lam = Array{Composites.laminate,1}(n_pt)
    precompinput = Array{PreComp.input,1}(n_pt)
    precompoutput = Array{PreComp.output,1}(n_pt)
    twist_d = twist*180.0/pi
    twistrate_d = PreComp.tw_rate(n_pt,sloc[1:n_pt],twist_d)
    # strainx = zeros(xafstrain)
    # strainy = zeros(yafstrain)
    # rot = zeros(2,2)
    for i = 1:n_pt
        # find leading edge
        lei = indmin(abs.(xaf[i,:]))
        # shift so leading edge is first
        xpc = circshift(xaf[i,:],-(lei-1))
        ypc = circshift(yaf[i,:],-(lei-1))

        # assemble input
        precompinput[i],mat[i],lam[i] = layup(nchord[i],twist_d[i],
        twistrate_d[i],xpc,ypc,lam_t[i,:],designparams,x_offset[i],orientation)

        # calculate composite properties: stiffness, mass, etc
        precompoutput[i] = PreComp.properties(precompinput[i])

        # rot[1,:] = [cos(-twist[i]),-sin(-twist[i])]
        # rot[2,:] = [sin(-twist[i]),cos(-twist[i])]
        #
        # xy = [xafstrain[i,:]*nchord[i],yafstrain[i,:]*nchord[i]]
        # xyrot = rot*xy
        # # determine locations at which to calculate strain
        # strainx[i,:] = xyrot[1] +xloc[i]
        # strainy[i,:] = xyrot[2] +zloc[i]
    end

    # Get ABD Matrices
    A = Array{Array{Float64,2},1}(n_pt)
    B = Array{Array{Float64,2},1}(n_pt)
    D = Array{Array{Float64,2},1}(n_pt)
    for i = 1:n_pt
        Q = Composites.getQ(mat[i])
        A[i],B[i],D[i] = Composites.getABD(lam[i],Q)
    end

    # Get buckling strains
    bucklingstrain = zeros(Float64,n_pt,2) #currently hardcoded for a single web
    _,bucklingstrain[:,1] = Composites.localbuckling(A,D,webloc*geom.normalchord)
    _,bucklingstrain[:,2] = Composites.localbuckling(A,D,(1-webloc)*geom.normalchord)

    # Assemble structural properties
     return compstructure(mat,lam,A,B,D,bucklingstrain,precompinput,precompoutput)

 end

function p_layup(normalchord,twist_d,twistrate_d,xaf,yaf,
    p_lam_t,designparams,p_x_offset,p_orientation)

    t1 = p_lam_t[1]#uni
    t2 = p_lam_t[2]#weave

    mat = []

    # Materials Input
    e1 = zeros(length(designparams.p_usedmaterials))
    e2 = zeros(length(designparams.p_usedmaterials))
    g12 = zeros(length(designparams.p_usedmaterials))
    anu12 = zeros(length(designparams.p_usedmaterials))
    density = zeros(length(designparams.p_usedmaterials))

    for i = 1:length(designparams.p_usedmaterials)
        matnames = designparams.plyprops.names
        idx = find(matnames -> matnames == designparams.p_usedmaterials[i],matnames)
        material = designparams.plyprops.plies[idx[1]]
        push!(mat,material)
        e1[i] = material.e1
        e2[i] = material.e2
        g12[i] = material.g12
        anu12[i] = material.nu12
        density[i] = material.rho
    end

    web1 = designparams.p_webloc[1]
    xsec_nodeU=[0.0;web1;1.0]

    n_laminaU=[2,2]
    n_pliesU=[1,1,
        1,1]

    t_lamU=[t1,t2,
        t1,t2]

    tht_lamU=[p_orientation[1],p_orientation[2],
        p_orientation[1],p_orientation[2]]
    mat_lamU=[1,2,
        1,2]

    # Lower surface
    xsec_nodeL = xsec_nodeU
    n_laminaL = n_laminaU
    n_pliesL = n_pliesU
    t_lamL = t_lamU
    tht_lamL = tht_lamU
    mat_lamL = mat_lamU

    # Web
    loc_web = [web1]
    n_laminaW = [2]
    n_pliesW = [1,1]
    t_lamW = [0.0,0.0]
    tht_lamW = [p_orientation[1],p_orientation[2]]
    mat_lamW = [1,2]

    leloc = p_x_offset

    # assemble output
    lam = Composites.laminate{Int64,Float64}(mat_lamL[1:div(end,2)],n_pliesL[1:div(end,2)],t_lamL[1:div(end,2)],tht_lamL[1:div(end,2)])

    precompinput = PreComp.input(
    normalchord,
    -twist_d,-twistrate_d,
    leloc,xaf,yaf,
    e1,e2,g12,anu12,density,
    xsec_nodeU,n_laminaU,n_pliesU,t_lamU,tht_lamU,mat_lamU,
    xsec_nodeL,n_laminaL,n_pliesL,t_lamL,tht_lamL,mat_lamL,
    loc_web,n_laminaW,n_pliesW,t_lamW,tht_lamW,mat_lamW)

    return precompinput,mat,lam
end

function w_layup(normalchord,twist_d,twistrate_d,xaf,yaf,
    w_lam_t,designparams,w_x_offset,w_orientation)

    t1 = w_lam_t[1]#uni
    t2 = w_lam_t[2]#weave
    t3 = w_lam_t[3]#wing foam
    t4 = w_lam_t[4]#web fabric
    t5 = w_lam_t[5]#web foam

    mat = []

    # Materials Input
    e1 = zeros(length(designparams.w_usedmaterials))
    e2 = zeros(length(designparams.w_usedmaterials))
    g12 = zeros(length(designparams.w_usedmaterials))
    anu12 = zeros(length(designparams.w_usedmaterials))
    density = zeros(length(designparams.w_usedmaterials))

    for i = 1:length(designparams.w_usedmaterials)
        matnames = designparams.plyprops.names
        idx = find(matnames -> matnames == designparams.w_usedmaterials[i],matnames)
        material = designparams.plyprops.plies[idx[1]]
        push!(mat,material)
        e1[i] = material.e1
        e2[i] = material.e2
        g12[i] = material.g12
        anu12[i] = material.nu12
        density[i] = material.rho
    end

    web1 = designparams.w_webloc[1]
    xsec_nodeU=[0.0;web1;1.0]

    n_laminaU=[5,5]
    n_pliesU=[1,1,1,1,1,
        1,1,1,1,1]

    t_lamU=[t1,t2,t3,t2,t1,
        t1,t2,t3,t2,t1]

    tht_lamU=[w_orientation[1],w_orientation[2],0.0,w_orientation[2],w_orientation[1],
        w_orientation[1],w_orientation[2],0.0,w_orientation[2],w_orientation[1]]
    mat_lamU=[1,2,3,2,1,
        1,2,3,2,1]

    # Lower surface
    xsec_nodeL = xsec_nodeU
    n_laminaL = n_laminaU
    n_pliesL = n_pliesU
    t_lamL = t_lamU
    tht_lamL = tht_lamU
    mat_lamL = mat_lamU

    # Web
    loc_web = [web1]
    n_laminaW = [3]
    n_pliesW = [1,1,1]
    t_lamW = [t4,t5,t4]
    tht_lamW = [w_orientation[3],0.0,w_orientation[3]]
    mat_lamW = [2,3,2]

    leloc = w_x_offset

    # assemble output
    lam = Composites.laminate{Int64,Float64}(mat_lamL[1:div(end,2)],n_pliesL[1:div(end,2)],t_lamL[1:div(end,2)],tht_lamL[1:div(end,2)])

    precompinput = PreComp.input(
    normalchord,
    -twist_d,-twistrate_d,
    leloc,xaf,yaf,
    e1,e2,g12,anu12,density,
    xsec_nodeU,n_laminaU,n_pliesU,t_lamU,tht_lamU,mat_lamU,
    xsec_nodeL,n_laminaL,n_pliesL,t_lamL,tht_lamL,mat_lamL,
    loc_web,n_laminaW,n_pliesW,t_lamW,tht_lamW,mat_lamW)

    return precompinput,mat,lam
end

function beam_strain_wrapper(Fp,n_sec,ctlparams,printiter,plots,xaf,yaf,xafstrain,yafstrain,twist,x_offset,aerocenter,nchord,structure,spanyloc,xloc,yloc,zloc)

    momarmx_scac = zeros(n_sec) #shear center (sc) to aero center (ac) moment arm
    momarmy_scac = zeros(n_sec)
    momarmx_cmtc = zeros(n_sec)
    momarmy_cmtc = zeros(n_sec)
    momarmx_sccm = zeros(n_sec)
    momarmy_sccm = zeros(n_sec)

    if (printiter%(ctlparams.printfreq)==0.0 && plots)
        figure("af")
        clf()
    end

    strainx = zeros(xafstrain) # for strain locations about the shear center
    strainy = zeros(yafstrain)
    strainx_vis = zeros(xafstrain)
    strainy_vis = zeros(yafstrain)
    rot = zeros(2,2)

    for i = 1:n_sec
        # Calculate aerodynamic center from precalculated x-aero center and chord line
        idxle = indmin(abs.(xaf[i,:]))
        idxte = indmax(abs.(xaf[i,:]))

        yle = yaf[i,idxle]
        yte = yaf[i,idxte]

        acx = aerocenter[i]
        acy = (yte-yle)/(1)*acx + yle #y/mx+b

        y_ac = (acx+x_offset[i])*nchord[i] #x and y are flipped in precomp
        x_ac = acy*nchord[i]
        x_sc = x_ac + structure.precompoutput[i].x_sc
        y_sc = y_ac + structure.precompoutput[i].y_sc
        x_tc = x_ac + structure.precompoutput[i].x_tc
        y_tc = y_ac + structure.precompoutput[i].y_tc
        x_cm = x_ac + structure.precompoutput[i].x_cm
        y_cm = y_ac + structure.precompoutput[i].y_cm

        #moment arms defined positive as in Fig. 13 in precomp manual with tc above cm, x and y are flipped as well in precomp
        momarmx_scac[i] = x_ac-y_sc
        momarmy_scac[i] = y_ac-x_sc

        momarmx_cmtc[i] = y_tc-y_cm
        momarmy_cmtc[i] = x_tc-x_cm

        momarmx_sccm[i] = y_cm-y_sc
        momarmy_sccm[i] = x_cm-x_sc

        #-------- Calculate Rotated Strain Locations about the Shear Center --------#
        rot[1,:] = [cos(-twist[i]),-sin(-twist[i])]
        rot[2,:] = [sin(-twist[i]),cos(-twist[i])]

        xy = [xafstrain[i,:]*nchord[i]+x_offset[i]*nchord[i]-y_sc,yafstrain[i,:]*nchord[i]-x_sc] #Offest by the shear center, precomp FOR is swapped
        xyrot = rot*xy
        # determine locations at which to calculate strain
        strainx[i,:] = xyrot[1]
        strainy[i,:] = xyrot[2]

        # add offset for visualization
        xy = [xafstrain[i,:]*nchord[i]-0.125*nchord[i],yafstrain[i,:]*nchord[i]] #Offest by the shear center, precomp FOR is swapped
        xyrot = rot*xy
        strainx_vis[i,:] = xyrot[1] + xloc[i]
        strainy_vis[i,:] = xyrot[2] + zloc[i]

        # Plots
        if (i%5==0 || i==n_sec) && (printiter%(ctlparams.printfreq)==0.0 && plots)

            figure("af$i")
            # clf()
            plot((xaf[i,:]+x_offset[i])*nchord[i],yaf[i,:]*nchord[i])
            if i==n_sec
                plot(y_ac,x_ac,"x",label = "Aero Center")
                plot(y_sc,x_sc,"D",label = "Shear Center")
                plot(y_tc,x_tc,".",label = "Mass Center")
                plot(y_cm,x_cm,"+",label = "Tension Center")
            else
                plot(y_ac,x_ac,"x")
                plot(y_sc,x_sc,"D")
                plot(y_tc,x_tc,".")
                plot(y_cm,x_cm,"+")
            end
            axis("equal")
            legend(loc="center left", bbox_to_anchor=(1, 0.5))
            pause(0.001)

            figure("strain_af$i")
            # clf()
            plot(strainx[i,:],strainy[i,:])
            plot(xafstrain[i,:]*nchord[i]+x_offset[i]*nchord[i],yafstrain[i,:]*nchord[i])
            if i==n_sec
                plot(y_ac-y_sc,x_ac-x_sc,"x",label = "Aero Center")
                plot(y_sc-y_sc,x_sc-x_sc,"D",label = "Shear Center")
                plot(y_tc-y_sc,x_tc-x_sc,".",label = "Mass Center")
                plot(y_cm-y_sc,x_cm-x_sc,"+",label = "Tension Center")

                plot(y_ac,x_ac,"x",label = "Aero Center")
                plot(y_sc,x_sc,"D",label = "Shear Center")
                plot(y_tc,x_tc,".",label = "Mass Center")
                plot(y_cm,x_cm,"+",label = "Tension Center")
            else
                plot(y_ac-y_sc,x_ac-x_sc,"x")
                plot(y_sc-y_sc,x_sc-x_sc,"D")
                plot(y_tc-y_sc,x_tc-x_sc,".")
                plot(y_cm-y_sc,x_cm-x_sc,"+")

                plot(y_ac,x_ac,"x")
                plot(y_sc,x_sc,"D")
                plot(y_tc,x_tc,".")
                plot(y_cm,x_cm,"+")
            end
            axis("equal")
            legend(loc="center left", bbox_to_anchor=(1, 0.5))
            pause(0.001)

        end

    end

    #-------- TRANSLATE CCBLADE LOADS TO PRECOMP/BEAM LOADS -------#
    nodes = length(spanyloc)
    elements = nodes - 1

    #Beam goes from root to tip
    Py = Fp[3,:] #Lift is z in aero, y in beam
    Pz = Fp[1,:] #Drag is x in aero, z in beam
    Px = Fp[2,:] #compression is y in aero, x in beam

    # --- extract the point forces from distributed in order to apply moments from ac to sc etc -----
    DOF = 6
    F = zeros(DOF*nodes)
    for i = 1:elements
        start = (i-1)*DOF  # (0, 0) start of matrix

        _, _, Fsub = BeamFEA.beam_matrix(spanyloc[i+1] - spanyloc[i], [0.0,0.0], [0.0,0.0], [0.0,0.0],
        [0.0,0.0], [0.0,0.0], [0.0,0.0], Px[i:i+1], Py[i:i+1], Pz[i:i+1])

        idx = start+1:start+2*DOF
        F[idx] += Fsub
    end

    # These are in the beam frame of reference
    EIy = zeros(nodes)
    EIz = zeros(nodes)
    EA = zeros(nodes)
    GJ = zeros(nodes)
    rhoA = zeros(nodes)
    rhoJ = zeros(nodes)
    # Px = zeros(nodes) # Already Calculated
    # Py = zeros(nodes)
    # Pz = zeros(nodes)
    Fx = zeros(nodes)
    Fy = zeros(nodes)
    Fz = zeros(nodes)
    Mx = zeros(nodes)
    My = zeros(nodes)
    Mz = zeros(nodes)
    kx = zeros(nodes)
    ky = zeros(nodes)
    kz = zeros(nodes)
    kthetax = zeros(nodes)
    kthetay = zeros(nodes)
    kthetaz = zeros(nodes)

    idxf = 1
    for i = 1:nodes
        # these are all distributed properties
        mass = structure.precompoutput[i].mass
        iyy = structure.precompoutput[i].flap_iner
        ixx = structure.precompoutput[i].lag_iner
        d_sccm2 =momarmx_sccm[i]^2 + momarmy_sccm[i]^2
        izz = ixx+iyy + mass*d_sccm2 #perpendicular axis theorem and parallel axis theorm

        EIy[i] = structure.precompoutput[i].ei_flap
        EIz[i] = structure.precompoutput[i].ei_lag
        EA[i] = structure.precompoutput[i].ea
        GJ[i] = structure.precompoutput[i].gj
        rhoA[i] = mass
        rhoJ[i] = izz
        # Px[i] = 0.0 # Already Calculated
        # Py[i] = 0.0
        # Pz[i] = 0.0
        Fx[i] = 0.0 #TODO:Extra Point Loads centripetal force here
        Fy[i] = 0.0
        Fz[i] = 0.0
        Mx[i] = F[idxf+3]*momarmy_scac[i] + F[idxf+5]*momarmx_scac[i]
        My[i] = 0.0 #TODO:centripetal here, but tension center not aligned with shear center so...?
        Mz[i] = 0.0 #TODO:centripetal here, but tension center not aligned with shear center so...?

        idxf += DOF
    end

    #fixed at hub (future work may allow axial twisting via a passive pitching device)
    kx[1] = Inf
    ky[1] = Inf
    kz[1] = Inf
    kthetax[1] = Inf
    kthetay[1] = Inf
    kthetaz[1] = Inf

    #-------- FEA ANALYSIS -------#
    delta, freq, V, K, M, F = BeamFEA.fea_analysis(spanyloc, EIy, EIz, EA, GJ, rhoA, rhoJ, Px, Py, Pz,
    Fx, Fy, Fz, Mx, My, Mz, kx, ky, kz, kthetax, kthetay, kthetaz)
    # println(minimum(delta))
    #TODO: use non-rigid strain calculation, i.e. composite curvature
    #TODO: include shear strain
    # get strain for each x strip along the blade
    strains = zeros(length(spanyloc),length(strainy[1,:]))
    for i = 1:length(strainy[1,:])
        strains[:,i], Nx, Vy, Vz, Tx, Myout, Mzout = BeamFEA.strain(spanyloc,
        strainy[:,i], strainx[:,i], EIy, EIz, EA,
        Px, Py, Pz, Fx, Fy, Fz, Mx, My, Mz, false)
    end

    mass_structure = trapz(spanyloc,rhoA)

    return strains, delta, freq, mass_structure,strainx_vis,strainy_vis
end

function bucklingcon(structure,strain,webloc,normalchord)
  # Used for constraining buckling
  sf = 1.5 #safety factor

  # buckling strain is negative, strain is always positive

  # NOT SMOOTH - will cause problems if webloc is design variable
  c_buckling = zeros(size(strain,1),size(strain,2))
  for i = 1:size(c_buckling,1)
    for iloc = 1:size(c_buckling,2)
      # if structure.strainlocx[i,iloc] < webloc*normalchord[i] #|| webloc == 0.0
          c_buckling[i,iloc] = (structure.bucklingstrain[i,1]+sf*strain[i,iloc])
      # else
      #     c_buckling[i,iloc] = (structure.bucklingstrain[i,2]+sf*strain[i,iloc])
      # end
    end
  end

  # # SMOOTH, but unnecessarily conservative
  # c_buckling = zeros(size(strain,1),size(strain,2),size(def.bucklingstrain,2)*2)
  # for i = 1:size(c_buckling,1)
  #   for iloc = 1:size(c_buckling,2)
  #     for ibuckle = 1:size(def.bucklingstrain,2)
  #       c_buckling[i,iloc,2*ibuckle-1] = (-def.bucklingstrain[i,ibuckle]-sf*strain[i,iloc])/def.bucklingstrain[i,ibuckle]
  #       c_buckling[i,iloc,2*ibuckle] = (-def.bucklingstrain[i,ibuckle]+sf*strain[i,iloc])/def.bucklingstrain[i,ibuckle]
  #     end
  #   end
  # end

  return c_buckling[:]/1E2
end

function stresscalc(strain::Array{Float64,2},shear::Array{Float64,1},
    mat::Composites.material,theta_d::Float64)
    # Get rotated material stiffness matrix
    qbar = Composites.getQ(mat,theta_d)
    # Determine Stresses from Strains in individual plys
    nnodes = size(strain,1)
    nloc = size(strain,2)
    plystress = zeros(nnodes,nloc,3)
    for i = 1:nnodes #Run entire length of wing #
        for iloc = 1:nloc # Run at every specified strain location
            # Convert shear stress to shear strain
            gam12 = shear[i]/qbar[3,3]-qbar[1,3]/qbar[3,3]*strain[i,iloc]
            # Calculate laminate stresses
            stress = qbar*[strain[i,iloc],0,gam12]
            # Transform stresses to ply axes
            plystress[i,iloc,:] = Composites.rotstress(stress,theta_d)
        end
    end
    return plystress
end #stress_calc

function stresscon(stress,mat::Composites.material,criteria = "tsaiwu";BVID=0.65,sf=1.5)
    # Used to constrain stress
    sigma1 = sf*stress[:,:,1]
    sigma2 = sf*stress[:,:,2]
    tau12 = sf*stress[:,:,3]
    sigma1tu = BVID*mat.xt
    sigma2tu = BVID*mat.yt
    sigma1cu = BVID*mat.xc
    sigma2cu = BVID*mat.yc
    tau12u = BVID*mat.s
    if criteria=="maxstress"
        c_stress = zeros(Float64,6,size(sigma1)...)
        for i = 1:size(sigma1,1)
            for j = 1:size(sigma1,2)
                c_stress[:,i,j],_ = Composites.maxstress(sigma1[i,j],sigma2[i,j],tau12[i,j],
                sigma1tu,sigma1cu,sigma2tu,sigma2cu,tau12u)
            end
        end
    elseif criteria=="tsaiwu"
        c_stress = zeros(Float64,size(sigma1)...)
        for i = 1:size(sigma1,1)
            for j = 1:size(sigma1,2)
                c_stress[i,j],_ = Composites.tsaiwu(sigma1[i,j],sigma2[i,j],tau12[i,j],
                sigma1tu,sigma1cu,sigma2tu,sigma2cu,tau12u)
                # maximum(c_stress)
            end
        end
    elseif criteria=="hashinrotem"
        c_stress = zeros(Float64,4,size(sigma1)...)
        for i = 1:size(sigma1,1)
            for j = 1:size(sigma1,2)
                c_stress[:,i,j],_ = Composites.hashinrotem(sigma1[i,j],sigma2[i,j],tau12[i,j],
                sigma1tu,sigma1cu,sigma2tu,sigma2cu,tau12u)
            end
        end
    else
        error("Specified failure criteria has no implementation")
    end
    return c_stress[:]-1.0
end #stresscon

function stress_wrapper(usedmaterials,plyprops,spanlocy,orientation,strain,shear)

    c_stress = []
    stress = [] #used in vtk output
    stresslayer = zeros(length(usedmaterials),length(spanlocy),length(strain[1,:]))
    j=1
    for i = 1:length(usedmaterials)

        if !contains(usedmaterials[i],"foam") #don't check stress on foam
            matnames = plyprops.names
            idx = find(matnames -> matnames == usedmaterials[i],matnames)

            material = plyprops.plies[idx[1]]
            orien = orientation[j]
            j+=1 #so we don't have to have extra design variables in orientation

            stressi = stresscalc(strain,shear,material,orien)

            if i==1 #only calc once since it's the same for the whole structure (assuming thin layups)
                stress = sqrt.(stressi[:,:,1].^2+stressi[:,:,2].^2+stressi[:,:,3].^2) #TODO break up or calculate von-mises etc
            end

            # determine contraint values
            c_stressi = stresscon(stressi,material)
            stresslayer[i,:,:] = c_stressi

            c_stress = [c_stress;c_stressi]
        end
    end

    return c_stress, stress, stresslayer
end #stress_wrapper

function spancoord(spanyloc,Rtip,dihedral,sweep,twist,x_offset,chord)
    # Define number of sections
    nsections = length(spanyloc)
    # Main  sections
    xloc = zeros(nsections)
    yloc = zeros(nsections)
    zloc = zeros(nsections)
    xloc[2:end] = cumsum((spanyloc[2:end]-spanyloc[1:(end-1)]).*cos.(dihedral[1:nsections-1]).*sin.(sweep[1:nsections-1]))
    yloc[2:end] = cumsum((spanyloc[2:end]-spanyloc[1:(end-1)]).*cos.(dihedral[1:nsections-1]).*cos.(sweep[1:nsections-1]))
    zloc[2:end] = cumsum((spanyloc[2:end]-spanyloc[1:(end-1)]).*sin.(dihedral[1:nsections-1]).*cos.(sweep[1:nsections-1]))
    # Find y length of  by stepping back spanylochrough let
    yend = Rtip
    # Scale appropriately to match span
    scaling = yend/yloc[nsections]
    xloc = xloc*scaling - cos.(twist).*x_offset.*chord
    yloc = yloc*scaling
    zloc = zloc*scaling + sin.(twist).*x_offset.*chord
    # Get structural spanwise parameter
    sloc = zeros(nsections)
    for i = 2:nsections
        sloc[i] = sloc[i-1]+sqrt((xloc[i]-xloc[i-1])^2.0+
        (yloc[i]-yloc[i-1])^2.0+(zloc[i]-zloc[i-1])^2.0)
    end

    # Assemble full  parameters
    sloc = sloc
    xloc = xloc
    yloc = yloc
    zloc = zloc

    return sloc,xloc,yloc,zloc
end

function chordlengths(chord,sweep)
    # Determine number of sections
    nsections = length(sweep) + 1
    # Normal Chord lengths
    normalchord = zeros(nsections)
    normalchord[1] = chord[1]*cos(sweep[1])
    for i = 2:(nsections-1)
        normalchord[i] = chord[i]*cos((sweep[i-1]+sweep[i])/2.0)
    end
    normalchord[end] = chord[end]*cos(sweep[end])
    return normalchord
end

function printmaxviol(c,name)
    if !isempty(find(c.>-1e-4))
        println(string(name),": ",maximum(c))
    end
end

#
# (private)
# trapezoidal integration
#
function trapz(x::Array{Float64,1}, y::Array{Float64,1})  # integrate y w.r.t. x
    integral = 0.0
    for i = 1:length(x)-1
        integral += (x[i+1]-x[i])*0.5*(y[i] + y[i+1])
    end
    return integral
end

function linear_interp(x_array::Array{Float64,1},y_array::Array{Float64,1},xnew::Float64)

    i = 1
    if xnew<minimum(x_array) # cap at max and min
        return y_array[indmin(x_array)]
    elseif xnew>maximum(x_array)
        return y_array[indmax(x_array)]
    else
        for i = 1:length(x_array)-1
            if x_array[i]<=xnew && x_array[i+1]>=xnew
                break
            end
        end

        fraction = (xnew - x_array[i])/(x_array[i+1] - x_array[i])
        return y_array[i] + fraction * (y_array[i+1] - y_array[i])
    end

end

function linear_interp(x_array::Array{Float64,1},y_array::Array{Float64,1},xnew::Array{Float64,1})
    ynew = zeros(xnew)
    for i = 1:length(xnew)
        ynew[i] = linear_interp(x_array,y_array,xnew[i])
    end

    return ynew
end


text:
function lazy_append(a::Grid,b::Grid)
  AppendedGrid(a,b)
end

struct AppendedGrid{Dc,Dp,A,B} <: Grid{Dc,Dp}
  a::A
  b::B
  function AppendedGrid( a::Grid{Dc,Dp}, b::Grid{Dc,Dp}) where {Dc,Dp}
    new{Dc,Dp,typeof(a),typeof(b)}(a,b)
  end
end

function get_node_coordinates(t::AppendedGrid)
  xa = get_node_coordinates(t.a)
  xb = get_node_coordinates(t.b)
  if xa === xb
    xa
  else
    lazy_append(xa,xb)
  end
end

function get_cell_node_ids(t::AppendedGrid)
  xa = get_node_coordinates(t.a)
  xb = get_node_coordinates(t.b)
  idsa = get_cell_node_ids(t.a)
  idsb = get_cell_node_ids(t.b)
  if xa === xb
    lazy_append(idsa,idsb)
  else
    nxa = Int32(length(xa))
    idsc = lazy_map(Broadcasting(i->i+nxa),idsb)
    lazy_append(idsa,idsc)
  end
end

function get_reffes(t::AppendedGrid)
  ra = get_reffes(t.a)
  rb = get_reffes(t.b)
  if ra === rb
    ra
  else
    vcat(ra,rb)
  end
end

function get_cell_type(t::AppendedGrid)
  ra = get_reffes(t.a)
  rb = get_reffes(t.b)
  a = get_cell_type(t.a)
  b = get_cell_type(t.b)
  if ra === rb
    lazy_append(a,b)
  else
    nra = Int8(length(ra))
    c = lazy_map(i->i+nra,b)
    lazy_append(a,c)
  end
end

function get_cell_coordinates(trian::AppendedGrid)
  a = get_cell_coordinates(trian.a)
  b = get_cell_coordinates(trian.b)
  lazy_append(a,b)
end

function get_cell_ref_coordinates(trian::AppendedGrid)
  a = get_cell_ref_coordinates(trian.a)
  b = get_cell_ref_coordinates(trian.b)
  lazy_append(a,b)
end

# In this case, we do not want a lazy_append since it will become difficult to
# compress / expand the reffes.
#function get_cell_reffe(trian::AppendedGrid)
#  a = get_cell_reffe(trian.a)
#  b = get_cell_reffe(trian.b)
#  lazy_append(a,b)
#end

function get_cell_shapefuns(trian::AppendedGrid)
  a = get_cell_shapefuns(trian.a)
  b = get_cell_shapefuns(trian.b)
  lazy_append(a,b)
end

function get_cell_map(trian::AppendedGrid)
  a = get_cell_map(trian.a)
  b = get_cell_map(trian.b)
  lazy_append(a,b)
end

function get_facet_normal(trian::AppendedGrid)
  cm = get_cell_map(trian)
  a = get_facet_normal(trian.a)
  b = get_facet_normal(trian.b)
  lazy_append(a,b)
end

function lazy_append(a::Triangulation,b::Triangulation)
  AppendedTriangulation(a,b)
end

struct AppendedTriangulation{Dc,Dp,A,B} <: Triangulation{Dc,Dp}
  a::A
  b::B
  function AppendedTriangulation(
    a::Triangulation{Dc,Dp}, b::Triangulation{Dc,Dp}) where {Dc,Dp}
    @assert get_background_model(a) === get_background_model(b)
    new{Dc,Dp,typeof(a),typeof(b)}(a,b)
  end
end

get_background_model(t::AppendedTriangulation) = get_background_model(t.a)

function get_grid(t::AppendedTriangulation)
  a = get_grid(t.a)
  b = get_grid(t.b)
  lazy_append(a,b)
end

function get_glue(t::AppendedTriangulation,::Val{D}) where D
  a = get_glue(t.a,Val(D))
  b = get_glue(t.b,Val(D))
  if a==nothing || b==nothing
    return nothing
  end
  lazy_append(a,b)
end

function lazy_append(a::FaceToFaceGlue,b::FaceToFaceGlue)
  tface_to_mface = lazy_append(a.tface_to_mface,b.tface_to_mface)
  tface_to_mface_map = lazy_append(a.tface_to_mface_map,b.tface_to_mface_map)
  mface_to_tface = nothing
  FaceToFaceGlue(tface_to_mface,tface_to_mface_map,mface_to_tface)
end

function get_cell_shapefuns(trian::AppendedTriangulation)
  a = get_cell_shapefuns(trian.a)
  b = get_cell_shapefuns(trian.b)
  lazy_append(a,b)
end

function get_cell_map(trian::AppendedTriangulation)
  a = get_cell_map(trian.a)
  b = get_cell_map(trian.b)
  lazy_append(a,b)
end

function get_facet_normal(trian::AppendedTriangulation)
  cm = get_cell_map(trian)
  a = get_facet_normal(trian.a)
  b = get_facet_normal(trian.b)
  lazy_append(a,b)
end

function move_contributions(
  cell_mat::AppendedArray, trian::AppendedTriangulation)

  if length(cell_mat.a) == num_cells(trian.a) && length(cell_mat.b) == num_cells(trian.b)
    a,ta = move_contributions(cell_mat.a,trian.a)
    b,tb = move_contributions(cell_mat.b,trian.b)
    return lazy_append(a,b), lazy_append(ta,tb)
  else
    return cell_mat, trian
  end
end

#function compress_contributions(cell_mat::AppendedArray,trian::AppendedTriangulation)
#  if length(cell_mat.a) == num_cells(trian.a) && length(cell_mat.b) == num_cells(trian.b)
#    a = compress_contributions(cell_mat.a,trian.a)
#    b = compress_contributions(cell_mat.b,trian.b)
#    return lazy_append(a,b)
#  else
#    return cell_mat
#  end
#end
#
#function compress_ids(cell_ids::AppendedArray,trian::AppendedTriangulation)
#  if length(cell_ids.a) == num_cells(trian.a) && length(cell_ids.b) == num_cells(trian.b)
#    a = compress_ids(cell_ids.a,trian.a)
#    b = compress_ids(cell_ids.b,trian.b)
#    return lazy_append(a,b)
#  else
#    return cell_ids
#  end
#end

ext: jl


text:
using Distributions, Turing, Random, Parameters
ProjDir = @__DIR__
cd(ProjDir)
include("LNR.jl")

@model model(data, x) = begin
    μ = Array{Real,1}(undef,Nr)
    β0 ~ Normal(0, 1)
    β1 ~ Normal(1, 1)
    μ = β0 .+ β1.*x
    s ~ Uniform(0, pi/2)
    σ = tan(s)
    for i in 1:length(data)
        data[i] ~ LNR(μ, σ, 0.0)
    end
end

Random.seed!(3423)
Nreps = 30
error_count = fill(0.0, Nreps)
Nr = 20
Nobs = 5

for i in 1:Nreps
    x = rand(Normal(0, 1), Nr)
    β0 = rand(Normal(0, .5))
    β1 = rand(Normal(1, .5))
    μ = β0 .+ β1.*x
    σ = rand(Uniform(.5, 2))
    dist = LNR(μ=μ, σ=σ, ϕ=0.0)
    data = rand(dist, Nobs)
    chain = sample(model(data, x), NUTS(1000, .8), 2000, discard_adapt = false, progress = false)
    error_count[i] = sum(get(chain,:numerical_error)[1])
end

describe(error_count)

max_stars_repo_name: itsdfish/Miscellaneous

max_issues_repo_name: itsdfish/Miscellaneous

max_forks_repo_name: itsdfish/Miscellaneous


text:
using FractionalDim
using LightXML

function read_xml_mesh(mesh::AbstractString)::Vector{Segment{3}}
    xroot = root(parse_file(mesh))
    mesh = first(child_elements(xroot))
    @assert attribute(mesh, "celltype") == "interval"

    dim = parse(attribute(mesh, "dim"))

    nodes, cells, _ = child_elements(mesh)

    nnodes = parse(attribute(nodes, "size"))
    X = zeros((dim, nnodes))

    get_coords = if dim == 2
        (e, col) -> begin
            X[dim*col + 1] = parse(attribute(e, "x"))
            X[dim*col + 2] = parse(attribute(e, "y"))
        end
    else
        (e, col) -> begin
            X[dim*col + 1] = parse(attribute(e, "x"))
            X[dim*col + 2] = parse(attribute(e, "y"))
            X[dim*col + 3] = parse(attribute(e, "z"))
        end
    end

    col = 0
    for node in child_elements(nodes)
        get_coords(node, col)
        col += 1
    end

    segments = Vector{Segment{dim}}()
    for cell in child_elements(cells)
        v0 = X[:, parse(attribute(cell, "v0")) + 1]
        v1 = X[:, parse(attribute(cell, "v1")) + 1]

        push!(segments, Segment{dim}(Point{dim}(v0), Point{dim}(v1)))
    end
    segments
end

# Read
maybe_fractal = read_xml_mesh("vasc_mesh.xml");

# And now estimate the dimension
counter = BoxCounter(maybe_fractal)
sizes, counts = Vector{Float64}(), Vector{Int}()

state = start(counter, 16)

rate, rate_lsq = Inf, Inf
while !done(counter, state)
    state = next(counter, state)
    size, count = last(state)
    push!(sizes, size)
    push!(counts, count)

    if length(sizes) > 1
        rate = -log(counts[end]/counts[end-1])/log(sizes[end]/sizes[end-1])
    end

    println((size, count, rate))
    writedlm("vasc_mesh_FD.txt", [sizes counts])
end

max_stars_repo_name: MiroK/fractional-dim


text:
# have to define several getindex methods to avoid ambiguities with other unit ranges
@propagate_inbounds function Base.getindex(a::AbstractIndex, i::AbstractUnitRange{<:Integer})
    return _getindex(a, to_index(a, i))
end
@propagate_inbounds function Base.getindex(a::AbstractIndex, i::Integer)
    return _getindex(a, to_index(a, i))
end
@propagate_inbounds function Base.getindex(a::AbstractIndex, i)
    return _getindex(a, to_index(a, i))
end

_getindex(idx::AbstractIndex, i::AbstractVector) = _maybe_index(idx, @inbounds(values(idx)[i]), i)
_getindex(idx::AbstractIndex, i) = @inbounds(values(idx)[i])

function _maybe_index(idx::Index{name}, vs::AbstractUnitRange{<:Integer}, i) where {name}
    return Index{name}(@inbounds(keys(idx)[i]), vs, UnkownUnique, LengthChecked)
end
function _maybe_index(idx::SimpleIndex{name}, vs::AbstractUnitRange{<:Integer}, i) where {name}
    return SimpleIndex{name}(vs)
end
# getindex of values promotes to a vector that can't be an index
_maybe_index(idx::Index{name}, vs::AbstractVector, i) where {name} = vs
_maybe_index(idx::SimpleIndex{name}, vs::AbstractVector, i) where {name} = vs


text:
"
 A Web socket server sample

 author: Atsushi Sakai
"

using HTTP

function main()
    println(PROGRAM_FILE, " start!!")

    HTTP.WebSockets.listen("127.0.0.1", UInt16(8081)) do ws
        while !eof(ws)
            data = readavailable(ws)
            write(ws, data)
            println("request received!!")

            for i in 1:10
                println("Say hey")
                sleep(3)
                write(ws, "Hey")
            end
        end
    end

    println(PROGRAM_FILE, " Done!!")
end

if abspath(PROGRAM_FILE) == @__FILE__
    @time main()
end

max_stars_repo_name: AtsushiSakai/JSONServerAndClient

max_issues_repo_name: AtsushiSakai/JSONServerAndClient

max_forks_repo_name: AtsushiSakai/JSONServerAndClient


text:
using Printf
using Dates: AbstractTime

maybe_int(t) = isinteger(t) ? Int(t) : t

"""
    prettytime(t)

Convert a floating point value `t` representing an amount of time in seconds to a more
human-friendly formatted string with three decimal places. Depending on the value of `t`
the string will be formatted to show `t` in nanoseconds (ns), microseconds (μs),
milliseconds (ms), seconds, minutes, hours, days, or years.
"""
function prettytime(t)
    # Modified from: https://github.com/JuliaCI/BenchmarkTools.jl/blob/master/src/trials.jl
    iszero(t) && return "0 seconds"

    t = maybe_int(t)

    if t < 1e-9
        # No point going to picoseconds, just return something readable.
        return @sprintf("%.3e seconds", t)
    elseif t < 1e-6
        value, units = t * 1e9, "ns"
    elseif t < 1e-3
        value, units = t * 1e6, "μs"
    elseif t < 1
        value, units = t * 1e3, "ms"
    elseif t < minute
        value = t
        units = value == 1 ? "second" : "seconds"
    elseif t < hour
        value = maybe_int(t / minute)
        units = value == 1 ? "minute" : "minutes"
    elseif t < day
        value = maybe_int(t / hour)
        units = value == 1 ? "hour" : "hours"
    elseif t < year
        value = maybe_int(t / day)
        units = value == 1 ? "day" : "days"
    else
        value = maybe_int(t / year)
        units = value == 1 ? "year" : "years"
    end

    if isinteger(value)
        return @sprintf("%d %s", value, units)
    else
        return @sprintf("%.3f %s", value, units)
    end
end

prettytime(dt::AbstractTime) = "$dt"


text:
module API

import tree_sitter_jll

using CEnum

const libtreesitter = tree_sitter_jll.libtreesitter_path

# Structs

const TREE_SITTER_LANGUAGE_VERSION = 11
const TREE_SITTER_MIN_COMPATIBLE_LANGUAGE_VERSION = 9
const TSSymbol = UInt16
const TSFieldId = UInt16

struct TSSymbolMetadata
    visible::Cint
    named::Cint
end

struct ANONYMOUS3_reduce
    symbol::TSSymbol
    dynamic_precedence::Int16
    child_count::UInt8
    production_id::UInt8
end

struct ANONYMOUS2_params
    reduce::ANONYMOUS3_reduce
end

@cenum TSParseActionType::UInt32 begin
    TSParseActionTypeShift = 0
    TSParseActionTypeReduce = 1
    TSParseActionTypeAccept = 2
    TSParseActionTypeRecover = 3
end

struct TSParseAction
    params::ANONYMOUS2_params
    type::TSParseActionType
end

struct TSParseActionEntry
    action::TSParseAction
end

struct TSLexMode
    lex_state::UInt16
    external_lex_state::UInt16
end

struct ANONYMOUS1_external_scanner
    states::Ptr{Cint}
    symbol_map::Ptr{TSSymbol}
    create::Ptr{Cvoid}
    destroy::Ptr{Cvoid}
    bool::Cvoid
    serialize::Ptr{Cvoid}
    deserialize::Ptr{Cvoid}
end

struct TSFieldMapSlice
    index::UInt16
    length::UInt16
end

struct TSFieldMapEntry
    field_id::TSFieldId
    child_index::UInt8
    inherited::Cint
end

struct TSLanguage
    version::UInt32
    symbol_count::UInt32
    alias_count::UInt32
    token_count::UInt32
    external_token_count::UInt32
    symbol_names::Ptr{Cstring}
    symbol_metadata::Ptr{TSSymbolMetadata}
    parse_table::Ptr{UInt16}
    parse_actions::Ptr{TSParseActionEntry}
    lex_modes::Ptr{TSLexMode}
    alias_sequences::Ptr{TSSymbol}
    max_alias_sequence_length::UInt16
    bool::Cvoid
    keyword_capture_token::TSSymbol
    external_scanner::ANONYMOUS1_external_scanner
    field_count::UInt32
    field_map_slices::Ptr{TSFieldMapSlice}
    field_map_entries::Ptr{TSFieldMapEntry}
    field_names::Ptr{Cstring}
    large_state_count::UInt32
    small_parse_table::Ptr{UInt16}
    small_parse_table_map::Ptr{UInt32}
    public_symbol_map::Ptr{TSSymbol}
end

const TSParser = Cvoid
const TSTree = Cvoid
const TSQuery = Cvoid
const TSQueryCursor = Cvoid

@cenum TSInputEncoding::UInt32 begin
    TSInputEncodingUTF8 = 0
    TSInputEncodingUTF16 = 1
end

@cenum TSSymbolType::UInt32 begin
    TSSymbolTypeRegular = 0
    TSSymbolTypeAnonymous = 1
    TSSymbolTypeAuxiliary = 2
end

struct TSPoint
    row::UInt32
    column::UInt32
end

struct TSRange
    start_point::TSPoint
    end_point::TSPoint
    start_byte::UInt32
    end_byte::UInt32
end

struct TSInput
    payload::Ptr{Cvoid}
    read::Ptr{Cvoid}
    encoding::TSInputEncoding
end

@cenum TSLogType::UInt32 begin
    TSLogTypeParse = 0
    TSLogTypeLex = 1
end

struct TSLogger
    payload::Ptr{Cvoid}
    log::Ptr{Cvoid}
end

struct TSInputEdit
    start_byte::UInt32
    old_end_byte::UInt32
    new_end_byte::UInt32
    start_point::TSPoint
    old_end_point::TSPoint
    new_end_point::TSPoint
end

struct TSNode
    context::NTuple{4, UInt32}
    id::Ptr{Cvoid}
    tree::Ptr{TSTree}
end

struct TSTreeCursor
    tree::Ptr{Cvoid}
    id::Ptr{Cvoid}
    context::NTuple{2, UInt32}
end

struct TSQueryCapture
    node::TSNode
    index::UInt32
end

struct TSQueryMatch
    id::UInt32
    pattern_index::UInt16
    capture_count::UInt16
    captures::Ptr{TSQueryCapture}
end

@cenum TSQueryPredicateStepType::UInt32 begin
    TSQueryPredicateStepTypeDone = 0
    TSQueryPredicateStepTypeCapture = 1
    TSQueryPredicateStepTypeString = 2
end

struct TSQueryPredicateStep
    type::TSQueryPredicateStepType
    value_id::UInt32
end

@cenum TSQueryError::UInt32 begin
    TSQueryErrorNone = 0
    TSQueryErrorSyntax = 1
    TSQueryErrorNodeType = 2
    TSQueryErrorField = 3
    TSQueryErrorCapture = 4
end

const ts_builtin_sym_end = 0
const TREE_SITTER_SERIALIZATION_BUFFER_SIZE = 1024

const TSStateId = UInt16

struct TSLexer
    lookahead::Int32
    result_symbol::TSSymbol
    advance::Ptr{Cvoid}
    mark_end::Ptr{Cvoid}
    get_column::Ptr{Cvoid}
    bool::Cvoid
end

struct ANONYMOUS4_external_scanner
    states::Ptr{Cint}
    symbol_map::Ptr{TSSymbol}
    create::Ptr{Cvoid}
    destroy::Ptr{Cvoid}
    bool::Cvoid
    serialize::Ptr{Cvoid}
    deserialize::Ptr{Cvoid}
end

# Functions

function ts_parser_new()
    ccall((:ts_parser_new, libtreesitter), Ptr{TSParser}, ())
end

function ts_parser_delete(self)
    ccall((:ts_parser_delete, libtreesitter), Cvoid, (Ptr{TSParser},), self)
end

function ts_parser_set_language(self, lang)
    ccall((:ts_parser_set_language, libtreesitter), Cint, (Ptr{TSParser}, Ptr{TSLanguage}), self, lang)
end

function ts_parser_language(self)
    ccall((:ts_parser_language, libtreesitter), Ptr{TSLanguage}, (Ptr{TSParser},), self)
end

function ts_parser_set_included_ranges()
    ccall((:ts_parser_set_included_ranges, libtreesitter), Cint, ())
end

function ts_parser_included_ranges(self, length)
    ccall((:ts_parser_included_ranges, libtreesitter), Ptr{TSRange}, (Ptr{TSParser}, Ptr{UInt32}), self, length)
end

function ts_parser_parse(self, old_tree, input)
    ccall((:ts_parser_parse, libtreesitter), Ptr{TSTree}, (Ptr{TSParser}, Ptr{TSTree}, TSInput), self, old_tree, input)
end

function ts_parser_parse_string(self, old_tree, string, length)
    ccall((:ts_parser_parse_string, libtreesitter), Ptr{TSTree}, (Ptr{TSParser}, Ptr{TSTree}, Cstring, UInt32), self, old_tree, string, length)
end

function ts_parser_parse_string_encoding(self, old_tree, string, length, encoding)
    ccall((:ts_parser_parse_string_encoding, libtreesitter), Ptr{TSTree}, (Ptr{TSParser}, Ptr{TSTree}, Cstring, UInt32, TSInputEncoding), self, old_tree, string, length, encoding)
end

function ts_parser_reset(self)
    ccall((:ts_parser_reset, libtreesitter), Cvoid, (Ptr{TSParser},), self)
end

function ts_parser_set_timeout_micros(self, timeout)
    ccall((:ts_parser_set_timeout_micros, libtreesitter), Cvoid, (Ptr{TSParser}, UInt64), self, timeout)
end

function ts_parser_timeout_micros(self)
    ccall((:ts_parser_timeout_micros, libtreesitter), UInt64, (Ptr{TSParser},), self)
end

function ts_parser_set_cancellation_flag(self, flag)
    ccall((:ts_parser_set_cancellation_flag, libtreesitter), Cvoid, (Ptr{TSParser}, Ptr{Cint}), self, flag)
end

function ts_parser_cancellation_flag()
    ccall((:ts_parser_cancellation_flag, libtreesitter), Ptr{Cint}, ())
end

function ts_parser_set_logger(self, logger)
    ccall((:ts_parser_set_logger, libtreesitter), Cvoid, (Ptr{TSParser}, TSLogger), self, logger)
end

function ts_parser_logger(self)
    ccall((:ts_parser_logger, libtreesitter), TSLogger, (Ptr{TSParser},), self)
end

function ts_parser_print_dot_graphs(self, file)
    ccall((:ts_parser_print_dot_graphs, libtreesitter), Cvoid, (Ptr{TSParser}, Cint), self, file)
end

function ts_tree_copy(self)
    ccall((:ts_tree_copy, libtreesitter), Ptr{TSTree}, (Ptr{TSTree},), self)
end

function ts_tree_delete(self)
    ccall((:ts_tree_delete, libtreesitter), Cvoid, (Ptr{TSTree},), self)
end

function ts_tree_root_node(self)
    ccall((:ts_tree_root_node, libtreesitter), TSNode, (Ptr{TSTree},), self)
end

function ts_tree_language(arg1)
    ccall((:ts_tree_language, libtreesitter), Ptr{TSLanguage}, (Ptr{TSTree},), arg1)
end

function ts_tree_edit(self, edit)
    ccall((:ts_tree_edit, libtreesitter), Cvoid, (Ptr{TSTree}, Ptr{TSInputEdit}), self, edit)
end

function ts_tree_get_changed_ranges(old_tree, new_tree, length)
    ccall((:ts_tree_get_changed_ranges, libtreesitter), Ptr{TSRange}, (Ptr{TSTree}, Ptr{TSTree}, Ptr{UInt32}), old_tree, new_tree, length)
end

# TODO: handle FILE.
# function ts_tree_print_dot_graph(arg1, arg2)
#     ccall((:ts_tree_print_dot_graph, libtreesitter), Cvoid, (Ptr{TSTree}, Ptr{FILE}), arg1, arg2)
# end

function ts_node_type(arg1)
    ccall((:ts_node_type, libtreesitter), Cstring, (TSNode,), arg1)
end

function ts_node_symbol(arg1)
    ccall((:ts_node_symbol, libtreesitter), TSSymbol, (TSNode,), arg1)
end

function ts_node_start_byte(arg1)
    ccall((:ts_node_start_byte, libtreesitter), UInt32, (TSNode,), arg1)
end

function ts_node_start_point(arg1)
    ccall((:ts_node_start_point, libtreesitter), TSPoint, (TSNode,), arg1)
end

function ts_node_end_byte(arg1)
    ccall((:ts_node_end_byte, libtreesitter), UInt32, (TSNode,), arg1)
end

function ts_node_end_point(arg1)
    ccall((:ts_node_end_point, libtreesitter), TSPoint, (TSNode,), arg1)
end

function ts_node_string(self)
    ccall((:ts_node_string, libtreesitter), Cstring, (TSNode,), self)
end

function ts_node_is_null(self)
    ccall((:ts_node_is_null, libtreesitter), Bool, (TSNode,), self)
end

function ts_node_is_named(self)
    ccall((:ts_node_is_named, libtreesitter), Bool, (TSNode,), self)
end

function ts_node_is_missing(self)
    ccall((:ts_node_is_missing, libtreesitter), Bool, (TSNode,), self)
end

function ts_node_is_extra(self)
    ccall((:ts_node_is_extra, libtreesitter), Bool, (TSNode,), self)
end

function ts_node_has_changes(self)
    ccall((:ts_node_has_changes, libtreesitter), Bool, (TSNode,), self)
end

function ts_node_has_error(self)
    ccall((:ts_node_has_error, libtreesitter), Bool, (TSNode,), self)
end

function ts_node_parent(self)
    ccall((:ts_node_parent, libtreesitter), TSNode, (TSNode,), self)
end

function ts_node_child(self, nth)
    ccall((:ts_node_child, libtreesitter), TSNode, (TSNode, UInt32), self, nth)
end

function ts_node_child_count(self)
    ccall((:ts_node_child_count, libtreesitter), UInt32, (TSNode,), self)
end

function ts_node_named_child(self, nth)
    ccall((:ts_node_named_child, libtreesitter), TSNode, (TSNode, UInt32), self, nth)
end

function ts_node_named_child_count(self)
    ccall((:ts_node_named_child_count, libtreesitter), UInt32, (TSNode,), self)
end

function ts_node_child_by_field_name(self, field_name, field_name_length)
    ccall((:ts_node_child_by_field_name, libtreesitter), TSNode, (TSNode, Cstring, UInt32), self, field_name, field_name_length)
end

function ts_node_child_by_field_id(self, field_id)
    ccall((:ts_node_child_by_field_id, libtreesitter), TSNode, (TSNode, TSFieldId), self, field_id)
end

function ts_node_next_sibling(self)
    ccall((:ts_node_next_sibling, libtreesitter), TSNode, (TSNode,), self)
end

function ts_node_prev_sibling(self)
    ccall((:ts_node_prev_sibling, libtreesitter), TSNode, (TSNode,), self)
end

function ts_node_next_named_sibling(self)
    ccall((:ts_node_next_named_sibling, libtreesitter), TSNode, (TSNode,), self)
end

function ts_node_prev_named_sibling(self)
    ccall((:ts_node_prev_named_sibling, libtreesitter), TSNode, (TSNode,), self)
end

function ts_node_first_child_for_byte(self, byte)
    ccall((:ts_node_first_child_for_byte, libtreesitter), TSNode, (TSNode, UInt32), self, byte)
end

function ts_node_first_named_child_for_byte(self, byte)
    ccall((:ts_node_first_named_child_for_byte, libtreesitter), TSNode, (TSNode, UInt32), self, byte)
end

function ts_node_descendant_for_byte_range(self, byte_from, byte_to)
    ccall((:ts_node_descendant_for_byte_range, libtreesitter), TSNode, (TSNode, UInt32, UInt32), self, byte_from, byte_to)
end

function ts_node_descendant_for_point_range(self, point_from, point_to)
    ccall((:ts_node_descendant_for_point_range, libtreesitter), TSNode, (TSNode, TSPoint, TSPoint), self, point_from, point_to)
end

function ts_node_named_descendant_for_byte_range(self, byte_from, byte_to)
    ccall((:ts_node_named_descendant_for_byte_range, libtreesitter), TSNode, (TSNode, UInt32, UInt32), self, byte_from, byte_to)
end

function ts_node_named_descendant_for_point_range(self, point_from, point_to)
    ccall((:ts_node_named_descendant_for_point_range, libtreesitter), TSNode, (TSNode, TSPoint, TSPoint), self, point_from, point_to)
end

function ts_node_edit(arg1, arg2)
    ccall((:ts_node_edit, libtreesitter), Cvoid, (Ptr{TSNode}, Ptr{TSInputEdit}), arg1, arg2)
end

function ts_node_eq(left, right)
    ccall((:ts_node_eq, libtreesitter), Bool, (TSNode, TSNode), left, right)
end

function ts_tree_cursor_new(node)
    ccall((:ts_tree_cursor_new, libtreesitter), TSTreeCursor, (TSNode,), node)
end

function ts_tree_cursor_delete(cursor)
    ccall((:ts_tree_cursor_delete, libtreesitter), Cvoid, (Ptr{TSTreeCursor},), cursor)
end

function ts_tree_cursor_reset(cursor, node)
    ccall((:ts_tree_cursor_reset, libtreesitter), Cvoid, (Ptr{TSTreeCursor}, TSNode), cursor, node)
end

function ts_tree_cursor_current_node(cursor)
    ccall((:ts_tree_cursor_current_node, libtreesitter), TSNode, (Ptr{TSTreeCursor},), cursor)
end

function ts_tree_cursor_current_field_name(cursor)
    ccall((:ts_tree_cursor_current_field_name, libtreesitter), Cstring, (Ptr{TSTreeCursor},), cursor)
end

function ts_tree_cursor_current_field_id(cursor)
    ccall((:ts_tree_cursor_current_field_id, libtreesitter), TSFieldId, (Ptr{TSTreeCursor},), cursor)
end

function ts_tree_cursor_goto_parent(cursor)
    ccall((:ts_tree_cursor_goto_parent, libtreesitter), Cint, (Ptr{TSTreeCursor},), cursor)
end

function ts_tree_cursor_goto_next_sibling(cursor)
    ccall((:ts_tree_cursor_goto_next_sibling, libtreesitter), Cint, (Ptr{TSTreeCursor},), cursor)
end

function ts_tree_cursor_goto_first_child(cursor)
    ccall((:ts_tree_cursor_goto_first_child, libtreesitter), Cint, (Ptr{TSTreeCursor},), cursor)
end

function ts_tree_cursor_goto_first_child_for_byte(cursor, byte)
    ccall((:ts_tree_cursor_goto_first_child_for_byte, libtreesitter), Int64, (Ptr{TSTreeCursor}, UInt32), cursor, byte)
end

function ts_tree_cursor_copy(cursor)
    ccall((:ts_tree_cursor_copy, libtreesitter), TSTreeCursor, (Ptr{TSTreeCursor},), cursor)
end

function ts_query_new(language, source, source_len, error_offset, error_type)
    ccall((:ts_query_new, libtreesitter), Ptr{TSQuery}, (Ptr{TSLanguage}, Cstring, UInt32, Ptr{UInt32}, Ptr{TSQueryError}), language, source, source_len, error_offset, error_type)
end

function ts_query_delete(query)
    ccall((:ts_query_delete, libtreesitter), Cvoid, (Ptr{TSQuery},), query)
end

function ts_query_pattern_count(query)
    ccall((:ts_query_pattern_count, libtreesitter), UInt32, (Ptr{TSQuery},), query)
end

function ts_query_capture_count(query)
    ccall((:ts_query_capture_count, libtreesitter), UInt32, (Ptr{TSQuery},), query)
end

function ts_query_string_count(query)
    ccall((:ts_query_string_count, libtreesitter), UInt32, (Ptr{TSQuery},), query)
end

function ts_query_start_byte_for_pattern(query, pattern)
    ccall((:ts_query_start_byte_for_pattern, libtreesitter), UInt32, (Ptr{TSQuery}, UInt32), query, pattern)
end

function ts_query_predicates_for_pattern(query, pattern_index, length)
    ccall((:ts_query_predicates_for_pattern, libtreesitter), Ptr{TSQueryPredicateStep}, (Ptr{TSQuery}, UInt32, Ptr{UInt32}), query, pattern_index, length)
end

function ts_query_capture_name_for_id(query, id, length)
    ccall((:ts_query_capture_name_for_id, libtreesitter), Cstring, (Ptr{TSQuery}, UInt32, Ptr{UInt32}), query, id, length)
end

function ts_query_string_value_for_id(query, id, length)
    ccall((:ts_query_string_value_for_id, libtreesitter), Cstring, (Ptr{TSQuery}, UInt32, Ptr{UInt32}), query, id, length)
end

function ts_query_disable_capture(arg1, arg2, arg3)
    ccall((:ts_query_disable_capture, libtreesitter), Cvoid, (Ptr{TSQuery}, Cstring, UInt32), arg1, arg2, arg3)
end

function ts_query_disable_pattern(arg1, arg2)
    ccall((:ts_query_disable_pattern, libtreesitter), Cvoid, (Ptr{TSQuery}, UInt32), arg1, arg2)
end

function ts_query_cursor_new()
    ccall((:ts_query_cursor_new, libtreesitter), Ptr{TSQueryCursor}, ())
end

function ts_query_cursor_delete(cursor)
    ccall((:ts_query_cursor_delete, libtreesitter), Cvoid, (Ptr{TSQueryCursor},), cursor)
end

function ts_query_cursor_exec(cursor, query, node)
    ccall((:ts_query_cursor_exec, libtreesitter), Cvoid, (Ptr{TSQueryCursor}, Ptr{TSQuery}, TSNode), cursor, query, node)
end

function ts_query_cursor_set_byte_range(arg1, arg2, arg3)
    ccall((:ts_query_cursor_set_byte_range, libtreesitter), Cvoid, (Ptr{TSQueryCursor}, UInt32, UInt32), arg1, arg2, arg3)
end

function ts_query_cursor_set_point_range(arg1, arg2, arg3)
    ccall((:ts_query_cursor_set_point_range, libtreesitter), Cvoid, (Ptr{TSQueryCursor}, TSPoint, TSPoint), arg1, arg2, arg3)
end

function ts_query_cursor_next_match(cursor, match)
    ccall((:ts_query_cursor_next_match, libtreesitter), Bool, (Ptr{TSQueryCursor}, Ptr{TSQueryMatch},), cursor, match)
end

function ts_query_cursor_remove_match(arg1, id)
    ccall((:ts_query_cursor_remove_match, libtreesitter), Cvoid, (Ptr{TSQueryCursor}, UInt32), arg1, id)
end

function ts_query_cursor_next_capture()
    ccall((:ts_query_cursor_next_capture, libtreesitter), Cint, ())
end

function ts_language_symbol_count(lang)
    ccall((:ts_language_symbol_count, libtreesitter), UInt32, (Ptr{TSLanguage},), lang)
end

function ts_language_symbol_name(lang, sym)
    ccall((:ts_language_symbol_name, libtreesitter), Cstring, (Ptr{TSLanguage}, TSSymbol), lang, sym)
end

function ts_language_symbol_for_name(lang, string, length, is_named)
    ccall((:ts_language_symbol_for_name, libtreesitter), TSSymbol, (Ptr{TSLanguage}, Cstring, UInt32, Cint), lang, string, length, is_named)
end

function ts_language_field_count(lang)
    ccall((:ts_language_field_count, libtreesitter), UInt32, (Ptr{TSLanguage},), lang)
end

function ts_language_field_name_for_id(lang, field_id)
    ccall((:ts_language_field_name_for_id, libtreesitter), Cstring, (Ptr{TSLanguage}, TSFieldId), lang, field_id)
end

function ts_language_field_id_for_name(lang, name, name_length)
    ccall((:ts_language_field_id_for_name, libtreesitter), TSFieldId, (Ptr{TSLanguage}, Cstring, UInt32), lang, name, name_length)
end

function ts_language_symbol_type(lang, sym)
    ccall((:ts_language_symbol_type, libtreesitter), TSSymbolType, (Ptr{TSLanguage}, TSSymbol), lang, sym)
end

# Language

import
    tree_sitter_bash_jll,
    tree_sitter_c_jll,
    tree_sitter_cpp_jll,
    tree_sitter_go_jll,
    tree_sitter_html_jll,
    tree_sitter_java_jll,
    tree_sitter_javascript_jll,
    tree_sitter_json_jll,
    tree_sitter_julia_jll,
    tree_sitter_php_jll,
    tree_sitter_python_jll,
    tree_sitter_ruby_jll,
    tree_sitter_rust_jll,
    tree_sitter_typescript_jll

const LANGUAGE_REGEX = r"^tree_sitter_(\w+)_jll$"
const LANGUAGES = Dict{Symbol,Tuple{Function,String}}()
for lang in filter(s -> occursin(LANGUAGE_REGEX, string(s)), names(@__MODULE__; imported=true))
    name = match(LANGUAGE_REGEX, string(lang))[1]
    func = Expr(:quote, Symbol(:tree_sitter_, name))
    lib = @eval $(Expr(:., lang, QuoteNode(Symbol(:libtreesitter_, name, :_path))))
    queries_dir = joinpath(getfield(@__MODULE__, lang).artifact_dir, "queries")
    f = @eval $(Symbol(:tree_sitter_, name))() = ccall(($func, $lib), Ptr{TSLanguage}, ())
    LANGUAGES[Symbol(name)] = (f, queries_dir)
end

function lang_ptr(n::Symbol)
    haskey(LANGUAGES, n) || error("unknown language '$n'")
    return get!(() -> LANGUAGES[n][1](), LANG_PTRS, n)
end
const LANG_PTRS = Dict{Symbol,Ptr{Cvoid}}()

function lang_queries(n::Symbol)
    haskey(LANGUAGES, n) || error("unknown language '$n'")
    _, dir = LANGUAGES[n]
    return get!(LANG_QUERIES, n) do
        dict = Dict{String,String}()
        # Custom/vendored queries that aren't available yet in the upstream
        # packages are loaded first. Queries for a particular language are
        # located in a subdirectory of src/queries named after the language.
        custom = joinpath(@__DIR__, "queries", string(n))
        dir = isdir(custom) ? custom : dir
        if isdir(dir)
            for file in readdir(dir)
                path = joinpath(dir, file)
                if isfile(path)
                    name, ext = splitext(file)
                    if ext == ".scm"
                        dict[name] = read(path, String)
                    end
                end
            end
        end
        return dict
    end
end
const LANG_QUERIES = Dict{Symbol,Dict{String,String}}()

end


text:
module WarpedImageSeries

#using Images, BlockRegistration, Interpolations
using Images, RegisterDeformation, Interpolations, MappedArrays

import Base: getindex, size, show
export WarpedSeriesView, warpedseries

struct WarpedSeriesView{TO,N,TI} <: AbstractArray{TO,N}
    img::AbstractArray{TI,N}
    tfms #can be anything that warpedview() accepts
    views #output of warpedview(img, tfm[i]) for each i
    shared_idxs #intersection of all sets of warped indices
end

#default to Float64 output precision
WarpedSeriesView(img::AbstractArray{T,N}, tfms, views, shared_idxs) where {T,N} =
    WarpedSeriesView{Float64,T,N}(img, tfms, views, shared_idxs)

#This is conservative in that we only include in the view pixels that are
#common to all images post-warping.
#Note: this index intersection may cause problems for general non-affine deformations;
#only tested on rigid deformations so far
function warpedseries(img::AbstractArray{T,N}, tfms, precision=Float64) where {T,N}
    @assert length(tfms) == size(img, N) #assume series is along last dimension
    sh_idxs = axes(img)
    vs = []
    colons = (fill(Colon(), N-1)...,)
    for i=1:length(tfms)
        v = warpedview(mappedarray(precision, view(img, colons..., i)), tfms[i])
        sh_idxs = ([intersect(x, y) for (x,y) in zip(sh_idxs, axes(v))]...,)
        push!(vs, v)
    end
    return WarpedSeriesView{precision,N,T}(img, tfms, vs, sh_idxs)
end

show(io::IO, A::T) where {T<:WarpedSeriesView} = print(io, "WarpedSeriesView of size $(size(A)), eltype $(eltype(A))\n")
show(io::IO, ::MIME"text/plain", A::T) where {T<:WarpedSeriesView} = show(io, A)

size(A::WarpedSeriesView{T,N}) where {T,N} = (map(length, A.shared_idxs)...,size(A.img, N))
getindex(A::WarpedSeriesView, inds...) = inner_index(A, last(inds), A.shared_idxs, Base.front(inds)...)
getindex(A::WarpedSeriesView, inds::CartesianIndex) = getindex(A, inds.I...) #inner_index(A, last(inds), A.shared_idxs, Base.front(inds)...)

inner_index(A::WarpedSeriesView, series_inds::T, inds_safe, inds_user...) where {T<:Union{Int, CartesianIndex}} = A.views[series_inds][calc_inds(inds_safe, inds_user...)...]

#when multiple stacks are included
function inner_index(A::WarpedSeriesView{T,N}, series_inds, inds_safe, inds_user...) where {T,N}
    if isa(series_inds, Colon)
        series_inds = last(axes(A.img))
    end
    prealloc = zeros(T, Base.front(_idx_shape(A, (inds_user...,2)))..., length(series_inds))
    inner_index!(prealloc, A, series_inds, inds_safe, inds_user...)
end

#Note: this is about half as efficient as it could be
function inner_index!(prealloc::AbstractArray{TA,N}, A::WarpedSeriesView{TA,N}, series_inds::T, inds_safe, inds_user...) where {TA, N, T<:AbstractUnitRange{Int}}
    colons = (fill(Colon(), N-1)...,) #performance hit?
    for (i, idx) in enumerate(series_inds)
        prealloc[colons...,i] =  A.views[idx][calc_inds(inds_safe, inds_user...)...]
    end
    return prealloc
end

calc_inds(inds_safe, inds_user...) = (first(inds_safe)[first(inds_user)], calc_inds(Base.tail(inds_safe), Base.tail(inds_user)...)...)
calc_inds(inds_safe::Tuple{}, inds_user::Tuple{}...) = ()
calc_inds(inds_safe, inds_user::Tuple{}...) = error("Too few indices passed to getindex")
calc_inds(inds_safe::Tuple{}, inds_user...) = begin @show inds_user; error("Too many indices passed to getindex") end

#note: duplicated from CachedSeries.jl
_idxs(A, dim, idxs) = idxs
_idxs(A, dim, idxs::Colon) = axes(A,dim)
function _idx_shape(A::AbstractArray, idxs)
    lens = zeros(Int, ndims(A))
    for i = 1:ndims(A)
        lens[i] = length(_idxs(A, i, idxs[i]))
    end
    return (lens...,)
end

end # module