text:
"""
    check_inputs(inputs, npcf)

Check the input array of particles, and raise an error if it is of the wrong format.
"""
function check_inputs(inputs, npcf::NPCF)
    if npcf.coords=="spherical"
        @assert (maximum(inputs[:,1])<2pi && minimum(inputs[:,1])>=0) "Phi coordinate must be in range [0, 2pi)"
        @assert (maximum(inputs[:,2])<pi && minimum(inputs[:,2])>=0) "Theta coordinate must be in range [0, pi)"
    else
        if npcf.periodic
            @assert maximum(inputs[:,1:npcf.D])<=npcf._boxsize "Maximum position must be less than boxsize"
            @assert minimum(inputs[:,1:npcf.D])>=0. "Minimum position must be lessgreater than zero"
        end
    end
    @assert length(inputs[1,:])==npcf.D+1 "Inputs should have (D+1) values per particles!"
end

"""
    sep(p12, npcf)

Compute the 1D separation of two particles with separation vector `p12` in Cartesian or Spherical geometries.
For spherical geometries, we return the cosine of the angular separation along the great circle connecting the two points.
In this case, we assume that the separation vector is the coordinates of the second point with the first at the origin.
"""
@inline function sep(p12,npcf::NPCF)
    if npcf.coords=="cartesian"
        return sqrt(sum(p12.^2))
    elseif npcf.coords=="spherical"
        return cos(p12[2])
    else
        return error("unrecognized coordinate system!")
    end
end

"""
    periodic_wrapping!(sep_array, npcf)

Add periodic wrapping to an array of separation vectors `sep_array`.
"""
@inline function periodic_wrapping!(sep_array,npcf::NPCF)
    for k in length(sep_array[:,1])
        for i in 1:npcf.D
            while sep_array[k,i]<-npcf._boxsize/2
                sep_array[k,i] += npcf._boxsize
            end
            while sep_array[k,i]>npcf._boxsize/2
                sep_array[k,i] -= npcf._boxsize
            end
        end
    end
    return sep_array
end

"""
    reorient_points(input, origin)

Rotate a set of input points (`input`) such that a given origin point (`origin`) is on the z-axis.
This applies only to points on the 2-sphere in spherical coordinates.
"""
@inline function reorient_points(input,origin)

    ## First generate a rotation matrix
    X=[sin(origin[2])cos(origin[1]),sin(origin[2])sin(origin[1]),cos(origin[2])]

    # Generate rotation axis
    axis=cross(X,[0,0,1])
    axis/=norm(axis)

    # Compute rotation matrix
    rotation_matrix = AngleAxis(origin[2],axis[1],axis[2],axis[3])

    ## Now rotate the input grid
    output = copy(input)
    for i in 1:length(input[:,1])
        phi = input[i,1]
        th = input[i,2]
        # Put into Cartesian coordinates
        cartesian_form = [sin(th)cos(phi),sin(th)sin(phi),cos(th)]
        # Apply rotation
        rot_cartesian = rotation_matrix*cartesian_form
        # Revert to polar coordinates
        rot_cartesian/=norm(rot_cartesian)
        output[i,1:2] = [atan(rot_cartesian[2],rot_cartesian[1]),acos(rot_cartesian[3])]
    end
    return output
end

"""
    return_bin(sep, npcf)

Compute the bin corresponding to a given separation, `sep`.
"""
@inline function return_bin(sep,npcf::NPCF)
    if sep<npcf.r_min
        return -1
    elseif sep>npcf.r_max
        return -1
    else
        return Int64(div(sep-npcf.r_min,npcf._dr)+1)
    end
end

"""
    computePlm(x, l, m)

Compute associated Legendre polynomial at `x` for integer or half-integer `l`,`m`
"""
function computePlm(x,l,m)
    eps = 1e-12
    1/gamma(1-m-eps)*_₂F₁(-l+eps,l+1+eps,1-m-eps,(1-x)/2)*((1+x)/(1-x))^(m/2)
end

"""
    computeYl(theta_array, l_array, NPCF)

Compute hyperspherical harmonic for input arrays of theta (`theta_array`) and ell (`l_array`)
"""
function computeYl(theta_array,l_array,npcf::NPCF)
    output = (-1)^(l_array[1])/sqrt(2*pi)*exp(im*l_array[1]*theta_array[1])
    for j in 2:npcf.D-1
        L = l_array[j]
        l = l_array[j-1]
        jPlL = sqrt((2L+j-1)/2*gamma(L+l+j-1)/gamma(L-l+1))*sin(theta_array[j])^((2-j)/2)*computePlm(cos(theta_array[j]),L+(j-2)/2,-(l+(j-2)/2))
        output *= jPlL
    end
    return output
end

"""
    hyperspherical_harmonic_single(r12, l_arr, npcf)

Return hyperspherical harmonics of particle, given position vector `r12` and single ell array `l_arr`.
For `D = 3` these are just the usual spherical harmonics.
In spherical coordinates, `r12` are the coordinates of the second particle, rotated such that the first particle is at the origin.
"""
function hyperspherical_harmonic_single(r12,l_arr,npcf::NPCF)

    if npcf.D==2
        # First define the angles
        if npcf.coords=="cartesian"
            theta1 = atan(r12[2],r12[1])
        elseif npcf.coords=="spherical"
            theta1 = r12[1]
        end
        return exp(im*l_arr[1]*theta1)/sqrt(2*pi)
    elseif npcf.D==3
        # First define the angles
        theta2 = acos(r12[3]) # i.e. theta
        theta1 = atan(r12[2],r12[1]) # i.e. phi

        # Compute P_l^m(theta) with normalization
        Plm = computePlmcostheta(theta2,l_arr[2],abs(l_arr[1]))[(l_arr[2],abs(l_arr[1]))]
        if l_arr[1]<0
            phase = (-1)^l_arr[1]
        else
            phase = 1
        end
        return phase/sqrt(2)*exp(im*l_arr[1]*theta1)*Plm
    elseif npcf.D==4
        # First define angles
        theta3 = acos(r12[4])
        theta2 = atan(sqrt(r12[2]^2+r12[1]^2),r12[3])
        theta1 = atan(r12[2],r12[1])

        # Assemble hyperspherical harmonics
        return computeYl([theta1,theta2,theta3],l_arr,npcf)
    else
        error("Not yet defined!")
    end
end

"""
    hyperspherical_harmonic(r12, npcf)

Return all hyperspherical harmonics of particle up to `lmax`, given Cartesian position vector `r12`.
For `D = 3`, these are just the usual spherical harmonics.
In spherical coordinates, `r12` are the coordinates of the second particle, rotated such that the first particle is at the origin.
"""
function hyperspherical_harmonic(r12,npcf::NPCF)

    if npcf.D==2
        # First define the angles
        if npcf.coords=="cartesian"
            theta1 = atan(r12[2],r12[1])
        elseif npcf.coords=="spherical"
            theta1 = r12[1]
        end

        # Assemble hyperspherical harmonics
        Ylm = zeros(ComplexF64,2*npcf.lmax+1)
        l_index = 1
        for l1 in -npcf.lmax:npcf.lmax
            Ylm[l_index] = exp(im*l1*theta1)/sqrt(2*pi)
            l_index += 1
        end
    elseif npcf.D==3
        # First define the angles
        theta2 = acos(r12[3]) # i.e. theta
        theta1 = atan(r12[2],r12[1]) # i.e. phi

        # Compute P_l^m(theta) with normalization
        Plm = computePlmcostheta(theta2,npcf.lmax)

        # Assemble hyperspherical harmonics
        Ylm = zeros(ComplexF64,(npcf.lmax+1)^2)
        l_index = 1
        for l_2 in 0:npcf.lmax # i.e. l
            for l_1 in -l_2:l_2 # i.e. m
                if l_1<0
                    phase = (-1)^l_1
                else
                    phase = 1
                end
                Ylm[l_index] = phase/sqrt(2)*exp(im*l_1*theta1)*Plm[(l_2,abs(l_1))]
                l_index += 1
            end
        end
    elseif npcf.D==4
        # First define angles
        theta3 = acos(r12[4])
        theta2 = atan(sqrt(r12[2]^2+r12[1]^2),r12[3])
        theta1 = atan(r12[2],r12[1])

        # Assemble hyperspherical harmonics
        Ylm = zeros(ComplexF64,Int64((npcf.lmax+1)*(npcf.lmax+2)*(2*npcf.lmax+3)/6))
        l_index = 1
        for l_3 in 0:npcf.lmax
            for l_2 in 0:l_3
                for l_1 in -l_2:l_2 # i.e. m
                    Ylm[l_index] = computeYl([theta1,theta2,theta3],[l_1,l_2,l_3],npcf)
                    l_index += 1
                end
            end
        end
    else
        error("Not yet defined!")
    end
    return Ylm
end

"""
    wigner9j(j1, j2, j3, j4, j5, j6, j7, j8, j9)

Define the Wigner 9j function in terms of the Wigner 6j functions
"""
function wigner9j(j1,j2,j3,j4,j5,j6,j7,j8,j9)
   output = 0
   for x in max(abs(j1-j9),abs(j4-j8),abs(j2-j6)):min(abs(j1+j9),abs(j4+j8),abs(j2+j6))
        tmp = (-1)^(2x)*(2x+1)*wigner6j(j1,j4,j7,j8,j9,x)*wigner6j(j2,j5,j8,j4,x,j6)*wigner6j(j3,j6,j9,x,j1,j2)
        output += tmp
    end
    return output
end

"""
    clebsch_gordan(l1_arr, l2_arr, L_arr, D)

Return the Clebsch-Gordan coefficient for three arrays of ells, `l1_arr`, `l2_arr` and `L_arr` in `D` dimensions.
"""
function clebsch_gordan(l1_arr,l2_arr,L_arr,D)
    if D==2
        (l1_arr+l2_arr)==L_arr ? 1. : 0.
    elseif D==3
        return (-1)^(-l1_arr[2]+l2_arr[2]-L_arr[1])*sqrt(2*L_arr[2]+1)*wigner3j(l1_arr[2],l2_arr[2],L_arr[2],l1_arr[1],l2_arr[1],-L_arr[1])
    elseif D==4
        pref = (-1)^(-l1_arr[2]+l2_arr[2]-L_arr[1])*(L_arr[3]+1)*sqrt((2*l1_arr[2]+1)*(2*l2_arr[2]+1)*(2*L_arr[2]+1))
        threej = wigner3j(l1_arr[2],l2_arr[2],L_arr[2],l1_arr[1],l2_arr[1],-L_arr[1])
        ninej = wigner9j(l1_arr[3]/2,l2_arr[3]/2,L_arr[3]/2,l1_arr[3]/2,l2_arr[3]/2,L_arr[3]/2,l1_arr[2],l2_arr[2],L_arr[2])
        return pref*threej*ninej
    else
        return error("not yet implemented!")
    end
end

"""
    basis_3pcf(l, r1, r2, npcf)

Return isotropic basis function of separation vectors `r1` and `r2` with primary index `l`.
For `D = 2`, this is a rescaled Chebyshev polynomial.
For `D = 3`, this is a rescaled Legendre polynomial.
For `D = 4`, we just use the spherical harmonic representation, for simplicity.
For spherical coordinates, we assume that `r1` and `r2` are rotated such that their polar angles give the opening angle of the triangle.
"""
function basis_3pcf(l,r1,r2,npcf::NPCF)
    if npcf.coords=="cartesian"
        mu = sum(r1.*r2)
    end

    if npcf.D==2
        if npcf.coords=="cartesian"
            return 1/(2*pi)*chebyshevt(l,mu)
        elseif npcf.coords=="spherical"
            phi = r1[1]-r2[1]
            return 1/(2*pi)*cos(l*phi)
        end
    elseif npcf.D==3
        return Pl(mu,l)*(-1)^l*sqrt(2*l+1)/(4*pi)
    elseif npcf.D==4
        if l==0
            return 1/(2pi^2)
        else
            # Define hyperspherical harmonics
            Yl1 = hyperspherical_harmonic(r1,npcf)
            Yl2 = hyperspherical_harmonic(r2,npcf)

            # Sum over intermediate states
            partial_sum = 0.
            for l_2 in 0:l
                for l_1 in -l_2:l_2
                    ax1 = Int64(1/6*l*(1+l)*(1+2l)+l_2*(l_2+1)+l_1+1)
                    ax2 = Int64(1/6*l*(1+l)*(1+2l)+l_2*(l_2+1)-l_1+1)
                    partial_sum += (-1)^(l_2-l_1)/(l+1)*real(Yl1[ax1]*Yl2[ax2])
                end
            end
            return partial_sum
        end
    else
        return error("not yet implemented!")
    end
end

"""
    basis_4pcf(l1, l2, l3, r1, r2, r3, npcf)

Return isotropic basis function of separation vectors `r1`, `r2` and `r3` with primary indices `l1`, `l2`, `l3`.
For `D = 2` it is simpler (and not inefficient) to work in polar coordinates, though a representation can be found in terms of dot products.
For `D = 3` and `D = 4`, we construct the basis from the spherical harmonics, which is straightforward, though more efficient representations can be found in terms of Cartesian dot products.
"""
function basis_4pcf(l1,l2,l3,r1,r2,r3,npcf::NPCF)

    if npcf.D==2
        # Cases with a zero are easy!
        if l1==0
            return (l2+l3==0)*basis_3pcf(abs(l2),r2,r3,npcf)/sqrt(2pi)
        elseif l2==0
            return (l1+l3==0)*basis_3pcf(abs(l1),r1,r3,npcf)/sqrt(2pi)
        elseif l3==0
            return (l1+l2==0)*basis_3pcf(abs(l1),r1,r2,npcf)/sqrt(2pi)
        end

        if npcf.coords=="cartesian"
            theta1 = atan(r1[2],r1[1])
            theta2 = atan(r2[2],r2[1])
            theta3 = atan(r3[2],r3[1])
        elseif npcf.coords=="spherical"
            theta1 = r1[1]
            theta2 = r2[1]
            theta3 = r3[1]
        end

        return cos(l1*theta1+l2*theta2+l3*theta3)/(2pi)^(3/2)
    elseif npcf.D==3
        # Cases with a zero are easy!
        if l1==0
            return (l2==l3)*basis_3pcf(l2,r2,r3,npcf)/sqrt(4pi)
        elseif l2==0
            return (l1==l3)*basis_3pcf(l1,r1,r3,npcf)/sqrt(4pi)
        elseif l3==0
            return (l1==l2)*basis_3pcf(l1,r1,r2,npcf)/sqrt(4pi)
        end

        output = 0.
        for l1_1 in -l1:l1
            Y1 = hyperspherical_harmonic_single(r1,[l1_1,l1],npcf)
            for l2_1 in -l2:l2
                Y2 = hyperspherical_harmonic_single(r2,[l2_1,l2],npcf)

                l3_1=-l1_1-l2_1
                if abs(l3_1)>l3; continue; end

                clebsch = clebsch_gordan([l1_1,l1],[l2_1,l2],[-l3_1,l3],npcf.D)
                clebsch *= (-1)^(l3-l3_1)/sqrt(2l3+1) # final CG coefficient
                if clebsch==0; continue; end

                Y3 = hyperspherical_harmonic_single(r3,[l3_1,l3],npcf)
                output += clebsch*real(Y1*Y2*Y3)
            end
        end
        return output

    elseif npcf.D==4
        # Cases with a zero are easy!
        if l1==0
            return (l2==l3)*basis_3pcf(l2,r2,r3,npcf)/(sqrt(2)pi)
        elseif l2==0
            return (l1==l3)*basis_3pcf(l1,r1,r3,npcf)/(sqrt(2)pi)
        elseif l3==0
            return (l1==l2)*basis_3pcf(l1,r1,r2,npcf)/(sqrt(2)pi)
        end

        output = 0.
        for l1_2 in 0:l1
            for l1_1  in -l1_2:l1_2
                Y1 = hyperspherical_harmonic_single(r1,[l1_1,l1_2,l1],npcf)
                for l2_2 in 0:l2
                    for l2_1 in -l2_2:l2_2
                        Y2 = hyperspherical_harmonic_single(r2,[l2_1,l2_2,l2],npcf)
                        for l3_2 in 0:l3
                            l3_1 = -l1_1-l2_1
                            if abs(l3_1)>l3_2; continue; end

                            clebsch = clebsch_gordan([l1_1,l1_2,l1],[l2_1,l2_2,l2],[-l3_1,l3_2,l3],npcf.D)
                            clebsch *= (-1)^(l3_2-l3_1)/(l3+1) # final CG coefficient
                            if clebsch==0.; continue; end

                            Y3 = hyperspherical_harmonic_single(r3,[l3_1,l3_2,l3],npcf)
                            output += clebsch*real(Y1*Y2*Y3)
                        end
                    end
                end
            end
        end
        return output
    end
end

"""
    basis_5pcf(l1, l2, l12, l3, l4, r1, r2, r3, r4, npcf)

Return isotropic basis function of `r1`, `r2`, `r3` and `r4` with primary indices `l1`, `l2`, `l12`, `l3`, `l4`.
For simplicity, we generate these from the isotropic basis functions, though representations can be found in terms of Cartesian dot products.
"""
function basis_5pcf(l1,l2,l12,l3,l4,r1,r2,r3,r4,npcf::NPCF)

    if npcf.D==2
        # Cases with a zero are easy!
        if l1==0
            return (l2+l3+l4==0)*basis_4pcf(l2,l3,l4,r2,r3,r4,npcf)/sqrt(2pi)
        elseif l2==0
            return (l1+l3+l4==0)*basis_4pcf(l1,l3,l4,r1,r3,r4,npcf)/sqrt(2pi)
        elseif l3==0
            return (l1+l2+l4==0)*basis_4pcf(l1,l2,l4,r1,r2,r4,npcf)/sqrt(2pi)
        elseif l4==0
            return (l1+l2+l3==0)*basis_4pcf(l1,l2,l3,r1,r2,r3,npcf)/sqrt(2pi)
        end

        if npcf.coords=="cartesian"
            theta1 = atan(r1[2],r1[1])
            theta2 = atan(r2[2],r2[1])
            theta3 = atan(r3[2],r3[1])
            theta4 = atan(r4[2],r4[1])
        elseif npcf.coords=="spherical"
            theta1 = r1[1]
            theta2 = r2[1]
            theta3 = r3[1]
            theta4 = r4[1]
        end
        return cos(l1*theta1+l2*theta2+l3*theta3+l4*theta4)/(2pi)^2
    elseif npcf.D==3
        # Cases with a zero are easy!
        if l1==0
            return (l2==l12)*basis_4pcf(l2,l3,l4,r2,r3,r4,npcf)/sqrt(4pi)
        elseif l2==0
            return (l1==l12)*basis_4pcf(l1,l3,l4,r1,r3,r4,npcf)/sqrt(4pi)
        elseif l3==0
            return (l4==l12)*basis_4pcf(l1,l2,l4,r1,r2,r4,npcf)/sqrt(4pi)
        elseif l4==0
            return (l3==l12)*basis_4pcf(l1,l2,l3,r1,r2,r3,npcf)/sqrt(4pi)
        end

        output = 0.
        for l1_1 in -l1:l1
            Y1 = hyperspherical_harmonic_single(r1,[l1_1,l1],npcf)
            for l2_1 in -l2:l2
                Y2 = hyperspherical_harmonic_single(r2,[l2_1,l2],npcf)
                l12_1 = l1_1+l2_1
                if abs(l12_1)>l12
                    continue
                end
                clebsch1 = clebsch_gordan([l1_1,l1],[l2_1,l2],[l12_1,l12],npcf.D)
                if clebsch1==0
                    continue
                end
                for l3_1 in -l3:l3
                    Y3 = hyperspherical_harmonic_single(r3,[l3_1,l3],npcf)
                    l4_1 = -l12_1-l3_1
                    if abs(l4_1)>l4
                        continue
                    end
                    Y4 = hyperspherical_harmonic_single(r4,[l4_1,l4],npcf)
                    clebsch2 = clebsch_gordan([l12_1,l12],[l3_1,l3],[-l4_1,l4],npcf.D)
                    if clebsch2==0
                        continue
                    end
                    clebsch3 = (-1)^(l4-l4_1)/sqrt(2l4+1) # this one is easy
                    output += clebsch1*clebsch2*clebsch3*real(Y1*Y2*Y3*Y4)
                end
            end
        end
        return output
    elseif npcf.D==4
        # Cases with a zero are easy!
        if l1==0
            return (l2==l12)*basis_4pcf(l2,l3,l4,r2,r3,r4,npcf)/(sqrt(2)pi)
        elseif l2==0
            return (l1==l12)*basis_4pcf(l1,l3,l4,r1,r3,r4,npcf)/(sqrt(2)pi)
        elseif l3==0
            return (l4==l12)*basis_4pcf(l1,l2,l4,r1,r2,r4,npcf)/(sqrt(2)pi)
        elseif l4==0
            return (l3==l12)*basis_4pcf(l1,l2,l3,r1,r2,r3,npcf)/(sqrt(2)pi)
        end

        output = 0.
        for l1_2 in 0:l1
            for l1_1  in -l1_2:l1_2
                Y1 = hyperspherical_harmonic_single(r1,[l1_1,l1_2,l1],npcf)
                for l2_2 in 0:l2
                    for l2_1 in -l2_2:l2_2
                        Y2 = hyperspherical_harmonic_single(r2,[l2_1,l2_2,l2],npcf)

                        for l12_2 in 0:l12
                            l12_1 = l1_1+l2_1
                            if abs(l12_1)>l12_2
                                continue
                            end

                            clebsch1 = clebsch_gordan([l1_1,l1_2,l1],[l2_1,l2_2,l2],[l12_1,l12_2,l12],npcf.D)
                            if clebsch1==0
                                continue
                            end

                            for l3_2 in 0:l3
                                for l3_1 in -l3_2:l3_2
                                    Y3 = hyperspherical_harmonic_single(r3,[l3_1,l3_2,l3],npcf)

                                    for l4_2 in 0:l4
                                        l4_1 = -l12_1-l3_1
                                        if abs(l4_1)>l4_2
                                            continue
                                        end
                                        Y4 = hyperspherical_harmonic_single(r4,[l4_1,l4_2,l4],npcf)
                                        clebsch2 = clebsch_gordan([l12_1,l12_2,l12],[l3_1,l3_2,l3],[-l4_1,l4_2,l4],npcf.D)
                                        if clebsch2==0
                                            continue
                                        end
                                        clebsch3 = (-1)^(l4_2-l4_1)/(l4+1) # this one is easy
                                        output += clebsch1*clebsch2*clebsch3*real(Y1*Y2*Y3*Y4)
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
        return output
    end
end

"""
    create_weight_array(npcf)

Create the array of angular weights for the NPCF calculation, using the Clebsch-Gordan coefficients.
"""
function create_weight_array(npcf::NPCF)
    if npcf.N==3
        if npcf.D==2
            weight_array = zeros(npcf.lmax+1)
            for l1_1 in 0:npcf.lmax
                # assuming l1_1 + l2_1 = 0, and l2_1<=0
                weight_array[l1_1+1] = 1
            end
        elseif npcf.D==3
            weight_array = zeros((npcf.lmax+1)^2)
            l_index = 1
            for l1_2 in 0:npcf.lmax
                for l1_1 in -l1_2:l1_2
                    # assuming l2_2 = l1_2 and l1_1 + l2_1 = 0
                    weight_array[l_index] = (-1)^(l1_2-l1_1)/sqrt(2l1_2+1)
                    l_index += 1
                end
            end
        elseif npcf.D==4
            weight_array = zeros(Int64((1+npcf.lmax)*(2+npcf.lmax)*(3+2npcf.lmax)/6))
            l_index = 1
            for l1_3 in 0:npcf.lmax
                for l1_2 in 0:l1_3
                    for l1_1 in -l1_2:l1_2
                        # assuming l1_3 = l2_3, l1_2 = l2_2 and l1_1 + l2_1 = 0
                        weight_array[l_index] = (-1)^(l1_2-l1_1)/(l1_3+1)
                        l_index += 1
                    end
                end
            end
        else
            return error("not yet implemented!")
        end
    elseif npcf.N==4
        if npcf.D==2
            weight_array = zeros(2npcf.lmax+1,2npcf.lmax+1)
            l1_ind = 1
            for l1_1 in -npcf.lmax:npcf.lmax
                l2_ind = 1
                for l2_1 in -npcf.lmax:npcf.lmax
                    if abs(l1_1+l2_1)>npcf.lmax || l1_1+l2_1<0
                        weight_array[l1_ind,l2_ind] = 0
                    else
                        weight_array[l1_ind,l2_ind] = 1
                    end
                    l2_ind += 1
                end
                l1_ind += 1
            end
        elseif npcf.D==3
            weight_array = zeros((npcf.lmax+1)^2,(npcf.lmax+1)^2,(npcf.lmax+1)^2)
            l1_ind = 1
            for l1_2 in 0:npcf.lmax
                for l1_1 in -l1_2:l1_2
                    l2_ind = 1
                    for l2_2 in 0:npcf.lmax
                        for l2_1 in -l2_2:l2_2
                            l3_ind = 1
                            for l3_2 in 0:npcf.lmax
                                for l3_1 in -l3_2:l3_2
                                    weight_array[l1_ind,l2_ind,l3_ind] = clebsch_gordan([l1_1,l1_2],[l2_1,l2_2],[-l3_1,l3_2],npcf.D)*clebsch_gordan([-l3_1,l3_2],[l3_1,l3_2],[0,0],npcf.D)
                                    l3_ind += 1
                                end
                            end
                            l2_ind += 1
                        end
                    end
                    l1_ind += 1
                end
            end
        elseif npcf.D==4
            ind_tot = Int64((npcf.lmax+1)*(npcf.lmax+2)*(2npcf.lmax+3)/6)
            weight_array = zeros(ind_tot,ind_tot,ind_tot)
            l1_ind = 1
            for l1_3 in 0:npcf.lmax
                for l1_2 in 0:l1_3
                    for l1_1 in -l1_2:l1_2
                        l2_ind = 1
                        for l2_3 in 0:npcf.lmax
                            for l2_2 in 0:l2_3
                                for l2_1 in -l2_2:l2_2
                                    l3_ind = 1
                                    for l3_3 in 0:npcf.lmax
                                        for l3_2 in 0:l3_3
                                            for l3_1 in -l3_2:l3_2
                                                weight_array[l1_ind,l2_ind,l3_ind] = clebsch_gordan([l1_1,l1_2,l1_3],[l2_1,l2_2,l2_3],[-l3_1,l3_2,l3_3],npcf.D)*clebsch_gordan([-l3_1,l3_2,l3_3],[l3_1,l3_2,l3_3],[0,0,0],npcf.D)
                                                l3_ind += 1
                                            end
                                        end
                                    end
                                    l2_ind += 1
                                end
                            end
                        end
                        l1_ind += 1
                    end
                end
            end
        end
    elseif npcf.N.==5
        if npcf.D==2
            weight_array = zeros(2npcf.lmax+1,2npcf.lmax+1,2npcf.lmax+1)
            l1_ind = 1
            for l1_1 in -npcf.lmax:npcf.lmax
                l2_ind = 1
                for l2_1 in -npcf.lmax:npcf.lmax
                    l3_ind = 1
                    for l3_1 in -npcf.lmax:npcf.lmax
                        if abs(l1_1+l2_1+l3_1)>npcf.lmax || l1_1+l2_1+l3_1<0
                            weight_array[l1_ind,l2_ind,l3_ind] = 0
                        else
                            weight_array[l1_ind,l2_ind,l3_ind] = 1
                        end
                        l3_ind += 1
                    end
                    l2_ind += 1
                end
                l1_ind += 1
            end
        elseif npcf.D==3
            weight_array = zeros((npcf.lmax+1)^2,(npcf.lmax+1)^2,2npcf.lmax+1,(npcf.lmax+1)^2,2npcf.lmax+1)
            l1_ind = 1
            for l1_2 in 0:npcf.lmax
                for l1_1 in -l1_2:l1_2
                    l2_ind = 1
                    for l2_2 in 0:npcf.lmax
                        for l2_1 in -l2_2:l2_2
                            l12_ind = 1
                            for l12_2 in 0:l1_2+l2_2
                                l12_1 = l1_1+l2_1
                                if abs(l12_1)>l12_2
                                    cb1 = 0
                                else
                                    cb1 = clebsch_gordan([l1_1,l1_2],[l2_1,l2_2],[l12_1,l12_2],npcf.D)
                                end
                                l3_ind = 1
                                for l3_2 in 0:npcf.lmax
                                    for l3_1 in -l3_2:l3_2
                                        l4_ind = 1
                                        for l4_2 in 0:npcf.lmax
                                            l4_1 = -l12_1-l3_1
                                            if abs(l4_1)>l4_2 || abs(l12_1)>l12_2
                                                cb2 = 0
                                            else
                                                cb2 = clebsch_gordan([l12_1,l12_2],[l3_1,l3_2],[-l4_1,l4_2],npcf.D)*clebsch_gordan([-l4_1,l4_2],[l4_1,l4_2],[0,0],npcf.D)
                                            end
                                            weight_array[l1_ind,l2_ind,l12_ind,l3_ind,l4_ind] = cb1*cb2
                                            l4_ind += 1
                                        end
                                        l3_ind += 1
                                    end
                                end
                                l12_ind += 1
                            end
                            l2_ind +=1
                        end
                    end
                    l1_ind +=1
                end
            end
        elseif npcf.D==4
            ind_tot = Int64((npcf.lmax+1)*(npcf.lmax+2)*(2npcf.lmax+3)/6) # l1, l2, l3 indices
            ind_tot2 = Int64((2npcf.lmax+1)*(2npcf.lmax+2)/2) # l12, l4 indices

            weight_array = zeros(ind_tot,ind_tot,ind_tot2,ind_tot,ind_tot2)
            l1_ind = 1
            for l1_3 in 0:npcf.lmax
                for l1_2 in 0:l1_3
                    for l1_1 in -l1_2:l1_2
                        l2_ind = 1
                        for l2_3 in 0:npcf.lmax
                            for l2_2 in 0:l2_3
                                for l2_1 in -l2_2:l2_2
                                    l12_ind = 1
                                    for l12_3 in 0:l1_3+l2_3
                                        for l12_2 in 0:l12_3
                                            l12_1 = l1_1+l2_1
                                            if abs(l12_1)>l12_2
                                                cb1 = 0
                                            else
                                                cb1 = clebsch_gordan([l1_1,l1_2,l1_3],[l2_1,l2_2,l2_3],[l12_1,l12_2,l12_3],npcf.D)
                                            end
                                            l3_ind = 1
                                            for l3_3 in 0:npcf.lmax
                                                for l3_2 in 0:l3_3
                                                    for l3_1 in -l3_2:l3_2
                                                        l4_ind = 1
                                                        for l4_3 in 0:npcf.lmax
                                                            for l4_2 in 0:l4_3
                                                                l4_1 = -l12_1-l3_1
                                                                if abs(l4_1)>l4_2 || abs(l12_1)>l12_2
                                                                    cb2 = 0
                                                                    cb3 = 0
                                                                else
                                                                    cb2 = clebsch_gordan([l12_1,l12_2,l12_3],[l3_1,l3_2,l3_3],[-l4_1,l4_2,l4_3],npcf.D)
                                                                    cb3 = clebsch_gordan([-l4_1,l4_2,l4_3],[l4_1,l4_2,l4_3],[0,0,0],npcf.D)
                                                                end
                                                                weight_array[l1_ind,l2_ind,l12_ind,l3_ind,l4_ind] = cb1*cb2*cb3
                                                                l4_ind += 1
                                                            end

                                                        end
                                                        l3_ind += 1
                                                    end
                                                end
                                            end
                                            l12_ind += 1
                                        end
                                    end
                                    l2_ind += 1
                                end
                            end
                        end
                        l1_ind += 1
                    end
                end
            end
        end
    end
    return weight_array
end

"""
    create_npcf_array(npcf)

Create the output array to store the NPCF calculation. The size of the array varies with the NPCF order and dimension.
"""
function create_npcf_array(npcf::NPCF)
    if npcf.N==2
        return zeros(Float64,npcf.nbins)
    elseif npcf.N==3
        return zeros(Float64,npcf.nbins,npcf.nbins,npcf.lmax+1)
    elseif npcf.N==4
        if npcf.D==2
            n_l = 0
            for l1 in -npcf.lmax:npcf.lmax
                for l2 in -npcf.lmax:npcf.lmax
                    if abs(l1+l2)>npcf.lmax || l1+l2<0
                        continue
                    end
                    n_l += 1
                end
            end
            return zeros(Float64,npcf.nbins,npcf.nbins,npcf.nbins,n_l)
        elseif npcf.D>=3
            n_l  = 0

            # Iterate over primary indices
            for l1 in 0:npcf.lmax
                for l2 in 0:npcf.lmax
                    for l3 in abs(l1-l2):min(l1+l2,npcf.lmax)
                        if (-1)^(l1+l2+l3)==1
                            # only even-parity NPCF contributions
                            n_l += 1
                        end
                    end
                end
            end
            return zeros(Float64,npcf.nbins,npcf.nbins,npcf.nbins,n_l)
        end
    elseif npcf.N==5
        if npcf.D==2
            n_l = 0
            for l1 in -npcf.lmax:npcf.lmax
                for l2 in -npcf.lmax:npcf.lmax
                    for l3 in -npcf.lmax:npcf.lmax
                        if abs(l1+l2+l3)>npcf.lmax || l1+l2+l3<0
                            continue
                        end
                        n_l += 1
                    end
                end
            end
            return zeros(Float64,npcf.nbins,npcf.nbins,npcf.nbins,npcf.nbins,n_l)
        elseif npcf.D>=3
            n_l  = 0

            # Iterate over primary indices
            for l1 in 0:npcf.lmax
                for l2 in 0:npcf.lmax
                    for l12 in abs(l1-l2):l1+l2
                        for l3 in 0:npcf.lmax
                            for l4 in abs(l12-l3):min(l12+l3,npcf.lmax)
                                if (-1)^(l1+l2+l3+l4)==1
                                    # only even-parity NPCF contributions
                                    n_l += 1
                                end
                            end
                        end
                    end
                end
            end
            return zeros(Float64,npcf.nbins,npcf.nbins,npcf.nbins,npcf.nbins,n_l)
        end
    else
        return error("not yet implemented!")
    end
end

"""
    accumulate_npcf_contribution(bins, Lambda, al, coupling_weights, npcf)

Accumulate the NPCF contribution for a single set of bins (`bins`) and primary ells (`Lambda`), given `al` and `coupling_weights` arrays.
This performs an outer product of (`N-1`) a_l arrays and a sum over intermediate ell (if `D>2`).
The bulk of the computation time is spent in this function for `N>3`.
"""
@inline function accumulate_npcf_contribution(bins,Lambda,al,coupling_weights,npcf::NPCF)
    contribution = 0

    if npcf.N==3
        b1,b2 = bins
        l1 = Lambda
        if npcf.D==2
            al1 = al[b1,npcf.lmax+l1+1]
            if al1==0; return 0.; end
            al2 = al[b2,npcf.lmax-l1+1]
            if al2==0; return 0.; end
            # taking real part and restricting to l1 >= 0
            contribution += coupling_weights[l1+1]*real(al1*al2)
        elseif npcf.D==3
            for l1_1 in -l1:l1
                ind1 = l1*(l1+1)+l1_1+1
                al1 = al[b1,ind1]
                if al1==0; continue; end
                ind2 = l1*(l1+1)-l1_1+1
                al2 = al[b2,ind2]
                if al2==0; continue; end
                contribution += coupling_weights[ind1]*real(al1*al2)
            end
        elseif npcf.D==4
            for l1_2 in 0:l1
                for l1_1 in -l1_2:l1_2
                    ind1 = Int64(1/6*l1*(1+l1)*(1+2l1)+l1_2*(l1_2+1)+l1_1+1)
                    al1 = al[b1,ind1]
                    if al1==0; continue; end
                    ind2 = Int64(1/6*l1*(1+l1)*(1+2l1)+l1_2*(l1_2+1)-l1_1+1)
                    al2 = al[b2,ind2]
                    if al2==0; continue; end
                    contribution += coupling_weights[ind1]*real(al1*al2)
                end
            end
        else
            return error("not yet implemented!")
        end

    elseif npcf.N==4
        b1,b2,b3 = bins
        if npcf.D==2
            l1_1, l2_1 = Lambda

            al1 = al[b1,l1_1+npcf.lmax+1]
            if al1==0; return 0.; end

            # Compute coupling matrix
            coupling = coupling_weights[l1_1+npcf.lmax+1,l2_1+npcf.lmax+1]
            if coupling==0; return 0.; end

            al2 = al[b2,l2_1+npcf.lmax+1]
            if al2==0; return 0.; end
            al3 = al[b3,-l1_1-l2_1+npcf.lmax+1]
            if al3==0; return 0.; end

            # Assemble contribution, taking real part and restricting to l1_1+l2_1>=0 implicitly
            contribution = real(al1*al2*al3)*coupling

        elseif npcf.D==3
            l1,l2,l3 = Lambda

            # Iterate over secondary ells (l3_1 is known)
            for l1_1 in -l1:l1
                al1 = al[b1,l1*(l1+1)+l1_1+1]
                if al1==0; continue; end
                for l2_1 in -l2:l2
                    al2 = al[b2,l2*(l2+1)+l2_1+1]
                    if al2==0; continue; end

                    l3_1 = -l1_1-l2_1
                    if abs(l3_1)>l3; continue; end
                    al3 = al[b3,l3*(l3+1)+l3_1+1]
                    if al3==0; continue; end

                    # Compute coupling matrix
                    coupling = coupling_weights[l1*(l1+1)+l1_1+1,l2*(l2+1)+l2_1+1,l3*(l3+1)+l3_1+1]
                    if coupling==0; continue; end

                    # Assemble contribution
                    contribution += real(al1*al2*al3)*coupling
                end
            end

        elseif npcf.D==4
            l1,l2,l3 = Lambda

            # Iterate over secondary ells (l3_1 is known)
            for l1_2 in 0:l1
                for l1_1 in -l1_2:l1_2
                    ind1 = Int64(1/6*l1*(1+l1)*(1+2l1)+l1_2*(l1_2+1)+l1_1+1)
                    al1 = al[b1,ind1]
                    if al1==0; continue; end

                    for l2_2 in 0:l2
                        for l2_1 in -l2_2:l2_2
                            ind2 = Int64(1/6*l2*(1+l2)*(1+2l2)+l2_2*(l2_2+1)+l2_1+1)
                            al2 = al[b2,ind2]
                            if al2==0; continue; end

                            for l3_2 in 0:l3
                                l3_1 = -l1_1-l2_1
                                if abs(l3_1)>l3_2; continue; end

                                ind3 = Int64(1/6*l3*(1+l3)*(1+2l3)+l3_2*(l3_2+1)+l3_1+1)

                                # Compute coupling matrix
                                coupling = coupling_weights[ind1,ind2,ind3]
                                if coupling==0; continue; end

                                al3 = al[b3,ind3]
                                if al3==0; continue; end

                                # Assemble contribution
                                contribution += real(al1*al2*al3)*coupling
                            end
                        end
                    end
                end
            end
        else
            return error("not yet implemented!")
        end

    elseif npcf.N==5
        b1,b2,b3,b4 = bins
        if npcf.D==2
            l1_1,l2_1,l3_1 = Lambda

            al1 = al[b1,l1_1+npcf.lmax+1]
            if al1==0; return 0.; end
            al2 = al[b2,l2_1+npcf.lmax+1]
            if al2==0; return 0.; end

            # Compute coupling matrix
            coupling = coupling_weights[l1_1+npcf.lmax+1,l2_1+npcf.lmax+1,l3_1+npcf.lmax+1]
            if coupling==0; return 0.; end

            al3 = al[b3,l3_1+npcf.lmax+1]
            if al3==0; return 0.; end
            al4 = al[b4,-l1_1-l2_1-l3_1+npcf.lmax+1]
            if al4==0; return 0.; end

            # Assemble contribution, taking real part and restricting to l1_1+l2_1>=0 implicitly
            contribution = real(al1*al2*al3*al4)*coupling

        elseif npcf.D==3
            l1,l2,l12,l3,l4=Lambda
            # Iterate over secondary ells (i.e. m) (l12_1, l4_1 are known)
            for l1_1 in -l1:l1
                al1 = al[b1,l1*(l1+1)+l1_1+1]
                if al1==0; continue; end
                for l2_1 in -l2:l2
                    al2 = al[b2,l2*(l2+1)+l2_1+1]
                    if al2==0; continue; end

                    l12_1 = l1_1+l2_1
                    if abs(l12_1)>l12; continue; end

                    for l3_1 in -l3:l3
                        al3 = al[b3,l3*(l3+1)+l3_1+1]
                        if al3==0; continue; end

                        l4_1 = -l12_1-l3_1
                        if abs(l4_1)>l4; continue; end

                        al4 = al[b4,l4*(l4+1)+l4_1+1]
                        if al4==0; continue; end

                        # Compute coupling matrix
                        coupling = coupling_weights[l1*(l1+1)+l1_1+1,l2*(l2+1)+l2_1+1,l12+1,l3*(l3+1)+l3_1+1,l4+1]
                        if coupling==0; continue; end

                        # Assemble contribution
                        contribution += real(al1*al2*al3*al4)*coupling
                    end
                end
            end

        elseif npcf.D==4
            l1,l2,l12,l3,l4=Lambda
            # Iterate over secondary ells (l12_1, l4_1 are known)
            for l1_2 in 0:l1
                for l1_1 in -l1_2:l1_2
                    ind1 = Int64(1/6*l1*(1+l1)*(1+2l1)+l1_2*(l1_2+1)+l1_1+1)
                    al1 = al[b1,ind1]
                    if al1==0; continue; end

                    for l2_2 in 0:l2
                        for l2_1 in -l2_2:l2_2
                            ind2 = Int64(1/6*l2*(1+l2)*(1+2l2)+l2_2*(l2_2+1)+l2_1+1)
                            al2 = al[b2,ind2]
                            if al2==0; continue; end

                            for l12_2 in 0:l12
                                l12_1 = l1_1 + l2_1
                                if abs(l12_1)>l12_2
                                    continue
                                end

                                for l3_2 in 0:l3
                                    for l3_1 in -l3_2:l3_2
                                        ind3 = Int64(1/6*l3*(1+l3)*(1+2l3)+l3_2*(l3_2+1)+l3_1+1)
                                        al3 = al[b3,ind3]
                                        if al3==0; continue; end

                                        for l4_2 in 0:l4
                                            l4_1 = -l12_1-l3_1
                                            if abs(l4_1)>l4_2; continue; end

                                            ind4 = Int64(1/6*l4*(1+l4)*(1+2l4)+l4_2*(l4_2+1)+l4_1+1)
                                            al4 = al[b4,ind4]
                                            if al4==0; continue; end

                                            # Compute coupling matrix
                                            coupling = coupling_weights[ind1,ind2,Int64(1/2*l12*(l12+1)+l12_2+1),ind3,Int64(1/2*l4*(l4+1)+l4_2+1)]
                                            if coupling==0; continue; end

                                            # Assemble contribution
                                            contribution += real(al1*al2*al3*al4)*coupling
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        else
            return error("not yet implemented!")
        end
    else
        return error("not yet implemented!")
    end
    return contribution
end;

"""
    normalize_npcf!(output, npcf)

Normalize the output NPCF measurements (`output`), using the binning parameters given in the `npcf` object.
This divides by the bin volume and the total volume of the space.
"""

function normalize_npcf!(output, npcf::NPCF)

    """Define volume of a single radial bin with index `b`."""
    function bin_vol(b)
        if npcf.D==2
            if npcf.coords=="spherical"
                return 2pi*npcf._dr
            else
                return pi*((npcf.r_min+npcf._dr*b)^2-(npcf.r_min+npcf._dr*(b-1))^2)
            end
        elseif npcf.D==3
            return 4pi/3*((npcf.r_min+npcf._dr*b)^3-(npcf.r_min+npcf._dr*(b-1))^3)
        elseif npcf.D==4
            return pi^2/2*((npcf.r_min+npcf._dr*b)^4-(npcf.r_min+npcf._dr*(b-1))^4)
        end
    end

    # Normalize by the bin volumes
    if npcf.N==2
        for b1 in 1:npcf.nbins
            normalization = bin_vol(b1)*npcf.volume
            output[b1]/=normalization
        end
    elseif npcf.N==3
        for b1 in 1:npcf.nbins
            for b2 in b1+1:npcf.nbins
                normalization = bin_vol(b1)bin_vol(b2)*npcf.volume
                output[b1,b2,:]/=normalization
            end
        end
    elseif npcf.N==4
        for b1 in 1:npcf.nbins
            for b2 in b1+1:npcf.nbins
                for b3 in b2+1:npcf.nbins
                    normalization = bin_vol(b1)bin_vol(b2)bin_vol(b3)*npcf.volume
                    output[b1,b2,b3,:]/=normalization
                end
            end
        end
    elseif npcf.N==5
        for b1 in 1:npcf.nbins
            for b2 in b1+1:npcf.nbins
                for b3 in b2+1:npcf.nbins
                    for b4 in b3+1:npcf.nbins
                        normalization = bin_vol(b1)bin_vol(b2)bin_vol(b3)bin_vol(b4)*npcf.volume
                        output[b1,b2,b3,b4,:]/=normalization
                    end
                end
            end
        end
    end
end;


text:
# Greetings! from path ~/.julia/config/startup.jl

# """
#     1. Generate tmp file: find $(pwd) -name \*.sac > ./tmp.txt
#     2. Run S0A_FFT.jl script
# """

############################################
## 1.Set parameters
############################################
# processes parameters
cores = 2                                               # total number of processes. cores number should be less than data-file number.
flag = true                                             # true of false; print computing time for debugging purpose.
precompile = false                                      # true or false; precompile or not.

# absolute path parameters
FFTDIR = expanduser("./FFT")                            # dir to store FFTData.

# select sac/mseed files or not
select_file = false                                     # true of false; whether specify sac/mseed files path.

# when select_file = ture, 'TMPTXT' works
FILETXT = expanduser("./file.txt")                      # sac/mseed files path.

# when select_file = false, 'DATADIR' 'select_sta' 'STATXT' and 'ncomp' work.
DATADIR = expanduser("./DATA2")                         # absolute dir where sac/mseed files are stored.
select_sta = "sta"                                      # "no" or "sta"; fft only works for selected net-sta base on STATXT.
STATXT = expanduser("./station_fft.txt")                # station infomation including: network-name and station-name.
channel_regular = "*CHN*"                               # channel(regular expression) of sac/mseed data; we will use 'glob(comp,DATADIR)' to select data in 'DATADIR' path.

# instrument parameters (The functionality has not been added!! please use [rm_resp = "no"])
rm_resp = "no"                                          # "no", "xml", "RESP", and "PZs"; select 'no' to not remove response and use "xml", "RESP", or "PZs" to remove response
respdir = expanduser("./XX.PZs")                        # path where resp files are located (required if rm_resp != "no")

# some control parameters
input_fmt = "sac"                                       # "sac" or "mseed".
time_norm = "one_bit"                                   # "no", "one_bit", "clip", "mute", and "smooth"; "no" for no normalization in time domain.
freq_norm = "whiten"                                    # "no", "whiten", "coherence", and "deconvolution"; "no" for no whitening, "whiten" for preserving only the phase component of the signal, "coherence" for smoothing absolute value of itself, and "deconvolution" for smoothing absolute value squared of itself.
cc_method = "CC"                                        # "CC" and "PCC"; FFT transformation depends on different cross-correlation methods, "CC" for pure cross correlation, and "PCC" for Phase cross-correlation, see [Ventosa et al., 2019] for "PCC".

# pre-processing parameters
freqmin = 0.5                                           # Float64; for 'bandpass!' of raw-data and "whiten" parameter
freqmax = 20.0                                          # Float64; for 'bandpass!' of raw-data and "whiten" parameter
cc_len = 60 * 2                                         # Int64; basic unit of data length for fft (sec)
cc_step = 60 * 1                                        # Int64; overlapping between each cc_len (sec)

# criteria for data selection in time_norm
factor_clip_std = 10.0                                  # Float64; 'clip' (default value), clipping the data whose std exceeds the 'factor_clip_std' to max/min mum value. || 中文: NoisePy与之对应的参数是 'max_over_std', 其默认值设置为10, 时间序列超过10倍均方差std的时间序列，把振幅设置成上下线（最大值），限制在一个范围内.
factor_mute = 3.0                                       # Float64; 'mute' (default value), mute the data whose amp exceed 'factor_mute' times the median to 0, for suppressing high-amplitude signals. || 中文: 时间序列超过3倍于波形的平均中值(使用 hibert 变换), 把振幅设置成0. median of envelope of `A` to find outliers.
time_half_win = 3                                       # Int64; 'smooth' (default value), the length(sample points) of the window for running mean normalization, using smooth(A, half_win) -> Smooth array A with half-window half_win (defaults to 3). || 中文: 将一段时间序列分成许多小的时窗，对每一个小的时窗内的波形进行滑动平均，NoisePy对应的参数是 "rma", 其窗口长度的默认值 smooth_N=10; 这也是bensen(2007)推荐的 time_norm 的方法.

# criteria for data selection in freq_norm
freq_half_win = 20                                      # Int64; 'coherence' and 'deconvolution' (default value), the length(sample points) of the window for running mean normalization, using coherence!(F,20) or deconvolution!(F,20) -> Smooth array F with half-window half_win (defaults to 20). || 中文: 频率域内的滑动平均，NoisePy对应的参数是 "rma", 其窗口长度的默认值 smooth_N=10; NoisePy 中 time_norm and freq_norm 使用的参数是同一个，都是smooth_N

############################################
# 2.run NoiseCC.FFT_CC
############################################
# a.using and import
# using Distributed
# NoiseCC_path = "/Users/yf/3.Project/2-1.Velocity_Change/SeisNoise_Code_YF/code_pmap_thread/"
# push!(LOAD_PATH, NoiseCC_path)
# @everywhere using NoiseCC # or @everywhere import NoiseCC

# b.include
using Distributed
addprocs(cores, topology = :master_worker)
@everywhere include("/Users/yf/3.Project/2-1.Velocity_Change/SeisNoise_Code_YF/code_pmap_thread/NoiseCC.jl")

NoiseCC.FFT(
    cores,
    flag,
    precompile,
    FFTDIR,
    select_file,
    FILETXT,
    DATADIR,
    select_sta,
    STATXT,
    input_fmt,
    time_norm,
    freq_norm,
    cc_method,
    channel_regular,
    rm_resp,
    respdir,
    freqmin,
    freqmax,
    cc_len,
    cc_step,
    factor_clip_std,
    factor_mute,
    time_half_win,
    freq_half_win,
)

ext: jl

max_stars_repo_name: OUCyf/NoiseCC

max_issues_repo_name: OUCyf/NoiseCC

max_forks_repo_name: OUCyf/NoiseCC


text:
"""
    MultiSampler

A sampler that makes it easy to sample from a vector of independent samplers.

##### Fields
- `samplers::T` : the independent samplers.
"""
struct MultiSampler{T<:AbstractVector{<:Sampleable}} <: Sampleable{Multivariate, Discrete}
    samplers::T
end

"""
    Base.show(io::IO, multisampler::MultiSampler) -> nothing

Show a `MultiSampler` in a human-friendly manner.
"""
function Base.show(io::IO, multisampler::MultiSampler)
    println(io, typeof(multisampler))
    println(io, " independent samplers: ")
    for sampler in multisampler.samplers
        println(io, "  ", sampler)
    end
    return nothing
end

"""
    MultiSampler(graph::AbstractGraph, nsamplers::Integer) -> MultiSampler

Construct and return a `MultiSampler` with `nsamplers` `NeighborSampler`s defined w.r.t.
`graph`.

##### Complexity
- Time complexity: `O(nsamplers)`.
"""
function MultiSampler(graph::AbstractGraph, nsamplers::Integer)
    return MultiSampler([NeighborSampler(graph, rand(1:nv(graph))) for _ in 1:nsamplers])
end

"""
    Base.rand(rng::AbstractRNG,
              multisampler::MultiSampler{<:AbstractVector{<:NeighborSampler}}) -> (
              Vector{Int})

Draw from each `NeighborSampler` independently and return a vector of the samples.

##### Complexity
- Time complexity: `O(nsamplers)`.

Use this for e.g., simulating multiple-agent random walks on a graph.
"""
function Base.rand(rng::AbstractRNG,
                   multisampler::MultiSampler{<:AbstractVector{<:NeighborSampler}})
    return [rand(sampler) for sampler in multisampler.samplers]
end

"""
    Base.rand(rng::AbstractRNG,
              multisampler::MultiSampler{<:AbstractVector{<:Sampleable}}) -> Vector{Float64}

Draw from each `sampler` independently and return a vector of the samples.

##### Complexity
- Time complexity: `O(nsamplers(multisampler))`.

Use this for e.g., simulating multiple-agent random walks on an integer lattice.
"""
function Base.rand(rng::AbstractRNG,
                   multisampler::MultiSampler{<:AbstractVector{<:Sampleable}})
    # convert to float for type stability
    return [float(rand(sampler)) for sampler in multisampler.samplers]
end

"""
    nsamplers(multisampler::MultiSampler)

Return the number of samplers.

##### Complexity
- Time complexity: `O(1)`.
"""
function nsamplers(multisampler::MultiSampler)
    return length(multisampler.samplers)
end

max_stars_repo_name: LoganDGraham/StochasticProcesses

max_issues_repo_name: LoganDGraham/StochasticProcesses

max_forks_repo_name: LoganDGraham/StochasticProcesses


text:
# DESCRIPTION
# Trims the aircraft dynamics at an arbitrary flight condition.
#
# INPUT
# - acfun: name of system dynamics function
# - x0: initial guess for trim state vector
# - u0: initial guess for trim control vector
# - targ_des: vector of trim targets
#
# OUTPUT
# - x0trim: trim state vector
# - u0trim: trim control vector
# - itrim: number of trim iterations
#
# ------------------------------------------------------------------------------

function trimmer(acfun,x0,u0,targ_des)
    #
    #x0trim=zeros(NSTATES)
    #u0trim=zeros(NCTRLS)
    #x0trim[:]=x0[:]
    #u0trim[:]=u0[:]
    x0trim=x0
    u0trim=u0
    #
    it=0
    conv=0
    trim_tol=5e-5
    itmax=80
    #
    @printf("trimming the aircraft\n")
    while ((it<itmax) && (conv==0))
        it=it+1
        # azimuthal average of xdot
        xdot0, y0=acfun(x0trim,u0trim)
        targvec=[xdot0;y0]
        targvec=targvec[TRIMTARG]
        targ_err=targvec-targ_des
        #
        @printf("%d          %f\n",it,norm(targ_err))
        #
        conv=1
        for k=1:NSTATES
            if (abs(targ_err[k])>TOL[k])
                conv=0
            end
        end
        #
        if (conv==0)
            #@show x0trim
            #@show u0trim
            A, B, C, D =linearize(acfun,x0trim,u0trim)
            #@show x0trim
            #@show u0trim
            #return
            # jacobian matrix
            Jac=[A B
                 C D]
            Jac2 = Jac[TRIMTARG,TRIMVARS]
            trimvec=[x0trim;u0trim]
            trimvec[TRIMVARS]=trimvec[TRIMVARS]-0.5*pinv(Jac2)*targ_err
            #trimvec[TRIMVARS]=trimvec[TRIMVARS]-Jac2\targ_err
            x0trim=trimvec[1:NSTATES]
            u0trim=trimvec[NSTATES+1:NSTATES+NCTRLS]
        end
    end
    #
    if (conv==0)
        @printf("warning: trim not achieved\n\n")
        itrim=0
    else
        @printf("successful trim\n\n")
        itrim=1
    end
    #
    return x0trim, u0trim, itrim
end

max_stars_repo_name: umbertosaetti/J-SimpleHel

max_issues_repo_name: umbertosaetti/J-SimpleHel

max_forks_repo_name: umbertosaetti/J-SimpleHel


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_leafletheatmap

"""
    ''_leafletheatmap(;kwargs...)

A LeafletHeatMap component.

Keyword arguments:
- `id` (String; optional)
- `blur` (Real; optional)
- `gradient` (Dict; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `max` (Real; optional)
- `minOpacity` (Real; optional)
- `points` (optional): . points has the following type: Array of lists containing elements 'lng', 'lat', 'weight'.
Those elements have the following types:
  - `lng` (Real; optional)
  - `lat` (Real; optional)
  - `weight` (Real; optional)s
- `radius` (Real; optional)
"""
function ''_leafletheatmap(; kwargs...)
        available_props = Symbol[:id, :blur, :gradient, :loading_state, :max, :minOpacity, :points, :radius]
        wild_props = Symbol[]
        return Component("''_leafletheatmap", "LeafletHeatMap", "feffery_leaflet_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: CNFeffery/feffery-leaflet-components

max_issues_repo_name: CNFeffery/feffery-leaflet-components

max_forks_repo_name: CNFeffery/feffery-leaflet-components


text:
abstract type NS2DProblem{T<:Real} end

findparams(::NS2DProblem{T}) where T = T

"""
    SimParams(Lx=2π, Ly=2π, nx=256, ny=256, ν=1e-4)

Parameters of the simulation. L denotes the size of the doubly
periodic computational domain, n the number of grid points to
discretize it in both direction and ν is the kinematic viscosity
of the working fluid.
"""
@with_kw struct SimParams
    # --> Dimension of the domain in the x-direction.
    Lx::Float64 = 2π ; @assert Lx > 0

    # --> Dimension of the domain in the y-direction.
    Ly::Float64 = 2π ; @assert Ly > 0

    # --> Number of grid points in the x-direction.
    nx::Int = 256 ; @assert mod(nx, 2) == 0

    # --> Number of grid points in the y-direction.
    ny::Int = 256 ; @assert mod(ny, 2) == 0

    # --> Viscosity.
    ν::Float64 = 1e-4 ; @assert ν > 0
end

#####################################
#####                           #####
#####     UNFORCED DYNAMICS     #####
#####                           #####
#####################################

mutable struct UnforcedProblem{T} <: NS2DProblem{T}
    # --> Number of grid points in the x-direction.
    nx::Int64

    # --> Number of grid points in the x-direction.
    ny::Int64

    # --> Viscosity.
    ν::T

    # --> Wavenumbers
    α::Array{T,1}
    β::Array{T,1}

    # --> Initial condition.
    ω₀::Array{Complex{T}, 2}

    # --> Working arrays for the computation of the advection term.
    adv::Array{Complex{T},2}
    ∂ω∂x::Array{Complex{T},2}
    ∂ω∂y::Array{Complex{T},2}
    u::Array{Complex{T},2}
    v::Array{Complex{T},2}

    # --> Working arrays for the computation of the dealiased advection term.
    advp::Array{Complex{T},2}
    ∂ω∂xp::Array{Complex{T},2}
    ∂ω∂yp::Array{Complex{T},2}
    up::Array{Complex{T},2}
    vp::Array{Complex{T},2}
end

function UnforcedProblem(p::SimParams, x::Array{Complex{T}, 2}) where T <: Real

    # --> Extract the parameters of the simulation.
    @unpack Lx, Ly, nx, ny, ν = p

    # --> Compute the wavenumbers.
    α = fftfreq(nx, nx/Lx) * 2π
    β = fftfreq(ny, ny/Ly) * 2π

    # --> Convert ν and α to whatever type T is.
    ν = convert(T, ν)
    α = convert(Array{T, 1}, α)
    β = convert(Array{T, 1}, β)

    # --> Template for the padded arrays.
    y = zeros(eltype(x), (3ny÷2, 3nx÷2))

    # --> Return the container.
    prob = UnforcedProblem(
        nx, ny, ν, α, β, x,
        zero(x), zero(x), zero(x), zero(x), zero(x),
        zero(y), zero(y), zero(y), zero(y), zero(y)
    )

    return prob
end

function unforced_dynamics!(dΩ::Array{Complex{T},2}, Ω::Array{Complex{T},2}, p::UnforcedProblem{T}, t::T) where T <: Real

    # --> Arrays for the computation of the advection term.
    adv, ∂ω∂x, ∂ω∂y, u, v = p.adv, p.∂ω∂x, p.∂ω∂y, p.u, p.v

    # --> Arrays for the dealiased computation.
    advp, ∂ω∂xp, ∂ω∂yp, up, vp = p.advp, p.∂ω∂xp, p.∂ω∂yp, p.up, p.vp

    # --> Compute all required quantities in spectral space.
    @inbounds for j = 1:p.nx
        @inbounds @simd for i = 1:p.ny
            # --> Misc.
            ω, α, β = Ω[i, j], p.α[j], p.β[i]
            k² = α^2 + β^2

            # --> Solve for the streamfunction.
            ψ = i == j == 1 ? zero(ω) : ω/k²

            # --> Diffusion term.
            dΩ[i, j] = -p.ν * k² * ω

            # --> Unpadded gradient of the vorticity field in spectral space.
            ∂ω∂x[i, j] = im * α * ω
            ∂ω∂y[i, j] = im * β * ω

            # --> Unpadded velocity in spectral space.
            u[i, j] = im * β * ψ
            v[i, j] = -im * α * ψ
        end
    end

    # # --> Enforce the Hermitian symmetry in the Fourier transform
    # #     since the loop above only acted on half the frequency space.
    # enforce_hermitian_symmetry!(dΩ)
    # enforce_hermitian_symmetry!(∂ω∂x), enforce_hermitian_symmetry!(∂ω∂y)
    # enforce_hermitian_symmetry!(u), enforce_hermitian_symmetry!(v)

    # --> Padding in spectral space.
    spectral_pad!(up, u), spectral_pad!(vp, v)
    spectral_pad!(∂ω∂xp, ∂ω∂x), spectral_pad!(∂ω∂yp, ∂ω∂y)

    # --> Compute the advection in physical space.
    ifft!(up), ifft!(vp), ifft!(∂ω∂xp), ifft!(∂ω∂yp)
    @. advp = up * ∂ω∂xp + vp * ∂ω∂yp

    # --> Chop the high-wavenumbers (anti-aliasing)
    spectral_chop!(adv, fft!(advp))

    # --> Compute the left-hand side.
    dΩ .-= adv
    dΩ[:, p.nx÷2+1] .= 0
    dΩ[p.ny÷2+1, :] .= 0

    return
end

function enforce_hermitian_symmetry!(x)

    n = size(x, 1)

    x[1, end:-1:n÷2+2] = conj(@view x[1, 2:n÷2])
    x[2:end, end:-1:n÷2+2] = conj(@view x[end:-1:2, 2:n÷2])

    return
end

function spectral_pad!(y, x)

    ny, nx = size(x)

    y[1:ny÷2, 1:nx÷2] = x[1:ny÷2, 1:nx÷2]
    y[1:ny÷2, end-nx÷2+1:end] = x[1:ny÷2, nx÷2+1:end]
    y[end-ny÷2+1:end, 1:nx÷2] = x[ny÷2+1:end, 1:nx÷2]
    y[end-ny÷2+1:end, end-nx÷2+1:end] = x[ny÷2+1:end, nx÷2+1:end]

    return
end

function spectral_chop!(y, x)

    ny, nx = size(y)

    y[1:ny÷2, 1:nx÷2] = x[1:ny÷2, 1:nx÷2]
    y[1:ny÷2, nx÷2+1:end] = x[1:ny÷2, end-nx÷2+1:end]
    y[ny÷2+1:end, 1:nx÷2] = x[end-ny÷2+1:end, 1:nx÷2]
    y[ny÷2+1:end, nx÷2+1:end] = x[end-ny÷2+1:end, end-nx÷2+1:end]

    return
end

rhs!(dΩ::Array{Complex{T},2}, Ω::Array{Complex{T},2}, p::UnforcedProblem{T}, t::T) where T <: Real = unforced_dynamics!(dΩ::Array{Complex{T},2}, Ω::Array{Complex{T},2}, p::UnforcedProblem{T}, t::T)

function simulate(nsprob::UnforcedProblem, T ; alg=Tsit5(), kwargs...)

    # --> Type of numbers used.
    N = findparams(nsprob)

    # --> Set the ODE Problem.
    tspan = (zero(N), convert(N, T))
    prob = ODEProblem(rhs!, copy(nsprob.ω₀), tspan, nsprob)

    # --> Integrate forward in time.
    sol = solve(prob, alg ; kwargs...)

    return sol
end


text:
"""
`RCropRatio <: ImageTransformation`

Description
============

Crops out the biggest possible area at some random position of
the given image, such that said sub-image satisfies the specified
aspect ratio (i.e. width divided by height).

Usage
======

    RCropRatio()

    RCropRatio(ratio = 1)

Arguments
==========

- **`ratio`** : The ratio of image height to image width that
the output image should satisfy

Methods
========

- **`transform`** : Applies the transformation to the given Image
or set of images

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)

Examples
========

    using Augmentor
    using TestImages

    # load an example image
    img = testimage("lena")

    # Randomly crops the image to a 2:1 aspect ratio
    img_new = transform(RCropRatio(2), img)

see also
=========

`ImageTransformation`, `transform`
"""
immutable RCropRatio <: ImageTransformation
    ratio::Float64

    function RCropRatio(ratio::Real)
        ratio > 0 || throw(ArgumentError("ratio has to be greater than 0"))
        new(Float64(ratio))
    end
end

RCropRatio(; ratio = 1.) = RCropRatio(ratio)

Base.showcompact(io::IO, op::RCropRatio) = print(io, "Crop to random $(op.ratio) aspect ratio.")
multiplier(::RCropRatio) = 1

function transform{T<:AbstractImage}(op::RCropRatio, img::T)
    w = width(img)
    h = height(img)

    nw = floor(Int, h * op.ratio)
    nh = floor(Int, w / op.ratio)
    nw = nw > 1 ? nw : 1
    nh = nh > 1 ? nh : 1

    result = if nw == w || nh == h
        img
    elseif nw < w
        x_max = w - nw + 1
        @assert x_max > 0
        x = rand(1:x_max)
        crop(img, x:(x+nw-1), 1:h)
    elseif nh < h
        y_max = h - nh + 1
        @assert y_max > 0
        y = rand(1:y_max)
        crop(img, 1:w, y:(y+nh-1))
    end::T
    _log_operation!(result, op)
end


text:
mutable struct UpdateInfo
	transform::Bool
	bound::Bool
end

mutable struct World <: BaseWorld
	parent::NodeObj
	children::Dict{Symbol, BaseObj}
	octree::Octree.Tree{NodeObj, Float32}
	updates::Dict{NodeObj, UpdateInfo}
	id::Symbol
	engine::Engine

	World(engine::Engine, aabb::AABB{Float32}, id::Symbol = :world) = new(NoNode(), Dict{Symbol, BaseObj}(), Octree.Tree(NodeObj, aabb), Dict{NodeObj, UpdateInfo}(), id, engine)
end

function init(world::World)
	add_event(world.engine, :render, world) do owner, event
		render(world)
	end
	add_event(world.engine, :update, world) do owner, event
		update_registered(world)
	end
	world
end

function init(engine::Engine, ::Type{World}, def::Dict{Symbol, Any})::World
	id = get_id!(def)
	bound = haskey(def, :bound) ? load_def(engine, def[:bound]) : AABB{Float32}([-256f0 256f0; -256f0 256f0; -256f0 256f0])
	world = init(World(engine, bound, id))
	init_children(engine, world, def)
	world
end

function done(world::World)
	remove_event(world, engine, :render)
end

function render_node(world::World, frustum::Convex{Float32}, node::Octree.Node{NodeObj}, nodeBound::AABB{Float32})
	outside(frustum, nodeBound) && return
	for obj in node.objects
		render(obj, world.engine)
	end
	for z = 1:2, y = 1:2, x = 1:2
		if !isa(node.subNodes[x, y, z], Octree.NullNode{NodeObj})
			subBound = Octree.getsubbound(nodeBound, [x, y, z])
			render_node(world, frustum, node.subNodes[x, y, z], subBound)
		end
	end
	nothing
end

function render(world::World)
	frustum = GRU.getfrustum(world.engine.renderer.camera)
	render_node(world, frustum, world.octree.root, world.octree.bound)
end

add_to_world(world::World, base::BaseObj) = nothing
function add_to_world(world::World, node::NodeObj)
	!haskey(node.children, :spatial) && return
	spatial = node.children[:spatial]
	worldMat = get_world_transform(node.parent)
	update_world_transform(spatial, worldMat)
	update_world_bound(spatial)
	add_node(world, node)
end

remove_from_world(world::World, base::BaseObj) = nothing
function remove_from_world(world::World, node::NodeObj)
	remove_node(world, node)
end

function add_node(world::World, node::NodeObj)
	bound = getbound(node)
	if !isempty(bound)
		Octree.add(world.octree, node, bound)
	end
end

function remove_node(world::World, node::NodeObj)
	bound = getbound(node)
	if !isempty(bound)
		Octree.remove(world.octree, node, bound)
	end
end

function register_transform_update(world::World, node::NodeObj)
	upd = get!(world.updates, node) do; UpdateInfo(false, false) end
	upd.transform = true
end

function register_bound_update(world::World, node::NodeObj)
	upd = get!(world.updates, node) do; UpdateInfo(false, false) end
	upd.bound = true
end

const defUpd = UpdateInfo(false, false)

update_transform(world::BaseWorld, leaf::BaseObj, worldMat::Matrix{Float32}) = nothing
function update_transform(world::World, node::NodeObj, worldMat::Matrix{Float32})
	!haskey(node.children, :spatial) && return
	remove_node(world, node)
	spatial = node.children[:spatial]
	upd = get(world.updates, node, defUpd)
	if upd.transform
		next_transform(spatial)
	end
	update_world_transform(spatial, worldMat)
	for child in node.children
		update_transform(world, child, get_world_transform(spatial))
	end
	if upd.bound
		next_bound(spatial)
		upd.bound = false
	end
	update_world_bound(spatial)
	add_node(world, node)
end

function update_registered(world::World)
	updateRoots = Set{NodeObj}()
	for (node, upd) in world.updates
		!upd.transform && continue
		root = node
		while true
			node = node.parent
			node == world && break
			if haskey(node.children, :spatial) && get(world.updates, node, defUpd).transform
				root = node
			end
		end
		push!(updateRoots, root)
	end
	for node in updateRoots
		update_transform(world, node, get_world_transform(node.parent))
	end
	for (node, upd) in world.updates
		!upd.bound && continue
		remove_node(world, node)
		spatial = node.children[:spatial]
		next_bound(spatial)
		update_world_bound(spatial)
		add_node(world, node)
	end
	empty!(world.updates)
end


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
Pkg.instantiate()

# ## Setup
using Dojo

# ## Mechanism
gravity = -9.81
timestep = 0.05
friction_coefficient = 0.8
damper = 1.0
env = get_environment(:quadruped,
    representation=:minimal,
    timestep=timestep,
    contact_body=true,
    gravity=gravity,
    friction_coefficient=friction_coefficient,
    damper=damper);

for j in env.mechanism.joints
    @show j.name
end

# ## Simulate
initialize!(env.mechanism, :quadruped,
    body_position=[0.0; 0.0; 0.5])

storage = simulate!(env.mechanism, 2.0,
    record=true,
    opts=SolverOptions(rtol=1.0e-4, btol=1.0e-4, verbose=false));

# ## Visualizer
open(env.vis)

# ## Visualize
visualize(env.mechanism, storage,
    vis=env.vis,
    show_contact=true)


text:
module TestMachines

using MLJBase
using Test
using Statistics
using ..Models

@load DecisionTreeRegressor

N=50
X = (a=rand(N), b=rand(N), c=rand(N));
y = 2*X.a - X.c + 0.05*rand(N);

train, test = partition(eachindex(y), 0.7);

tree = DecisionTreeRegressor(max_depth=5)

t = machine(tree, X, y)
@test_logs (:info, r"Training") fit!(t)
@test_logs (:info, r"Training") fit!(t, rows=train)
@test_logs (:info, r"Not retraining") fit!(t, rows=train)
@test_logs (:info, r"Training") fit!(t)
t.model.max_depth = 1
# MLJBase.recursive_setproperty!(t, :(model.max_depth),  1)
@test_logs (:info, r"Updating") fit!(t)

predict(t, selectrows(X,test));
@test rms(predict(t, selectrows(X, test)), y[test]) < std(y)

mach = machine(ConstantRegressor(), X, y)
@test_logs (:info, r"Training") fit!(mach)
yhat = predict_mean(mach, X);

n = nrows(X)
@test rms(yhat, y) ≈ std(y)*sqrt(1 - 1/n)

# test an unsupervised univariate case:
mach = machine(UnivariateStandardizer(), float.(1:5))
@test_logs (:info, r"Training") fit!(mach)
@test isempty(params(mach))

# test a frozen Machine
stand = machine(Standardizer(), source((x1=rand(10),)))
freeze!(stand)
@test_logs (:warn, r"not trained as it is frozen\.$") fit!(stand)

@testset "warnings" begin
    @test_throws DimensionMismatch machine(tree, X, y[1:end-1])
    @test_logs((:warn, r"The scitype of `y`"),
               machine(tree, X, categorical(1:N)))
    @test_logs((:warn, r"The scitype of `X`"),
               machine(tree, (x=categorical(1:N),), y))
end

@testset "weights" begin
    yraw = ["Perry", "Antonia", "Perry", "Skater"]
    X = (x=rand(4),)
    y = categorical(yraw)
    w = [2, 3, 2, 5]

    # without weights:
    mach = machine(ConstantClassifier(), X, y)
    fit!(mach)
    d1 = predict(mach, X)[1]
    d2 = MLJBase.UnivariateFinite([y[1], y[2], y[4]], [0.5, 0.25, 0.25])
    @test all([pdf(d1, c) ≈ pdf(d2, c) for c in MLJBase.classes(d1)])

    # with weights:
    mach = machine(ConstantClassifier(), X, y, w)
    fit!(mach)
    d1 = predict(mach, X)[1]
    d2 = MLJBase.UnivariateFinite([y[1], y[2], y[4]], [1/3, 1/4, 5/12])
    @test all([pdf(d1, c) ≈ pdf(d2, c) for c in MLJBase.classes(d1)])
end

mutable struct Scale <: MLJBase.Static
    scaling::Float64
end

function MLJBase.transform(s::Scale, _, X)
    X isa AbstractVecOrMat && return X * s.scaling
    MLJBase.table(s.scaling * MLJBase.matrix(X), prototype=X)
end

function MLJBase.inverse_transform(s::Scale, _, X)
    X isa AbstractVecOrMat && return X / s.scaling
    MLJBase.table(MLJBase.matrix(X) / s.scaling, prototype=X)
end

@testset "static transformer machines" begin
    s = Scale(2)
    X = randn(2, 3)
    Xt = MLJBase.table(X)

    @test_throws ArgumentError machine(Scale(2), X)
    @test_throws ArgumentError machine(Scale(2), source(X))

    mach = machine(Scale(2))
    fit!(mach) # no-op
    R  = transform(mach, X)
    IR = inverse_transform(mach, R)
    @test IR ≈ X
end

@testset "serialization" begin
    model = @load DecisionTreeRegressor

    X = (a = Float64[98, 53, 93, 67, 90, 68],
         b = Float64[64, 43, 66, 47, 16, 66],)
    Xnew = (a = Float64[82, 49, 16],
            b = Float64[36, 13, 36],)
    y =  [59.1, 28.6, 96.6, 83.3, 59.1, 48.0]

    mach =machine(model, X, y)
    filename = joinpath(@__DIR__, "machine.jlso")
    io = IOBuffer()
    @test_throws Exception MLJBase.save(io, mach; compression=:none)

    fit!(mach)
    report = mach.report
    pred = predict(mach, Xnew)
    MLJBase.save(io, mach; compression=:none)
    # commented out for travis testing:
    # MLJBase.save(filename, mach)

    # test restoring data from filename:
    m = machine(filename)
    p = predict(m, Xnew)
    @test m.model == model
    @test m.report == report
    @test p ≈ pred
    m = machine(filename, X, y)
    fit!(m)
    p = predict(m, Xnew)
    @test p ≈ pred

    # test restoring data from io:
    seekstart(io)
    m = machine(io)
    p = predict(m, Xnew)
    @test m.model == model
    @test m.report == report
    @test p ≈ pred
    seekstart(io)
    m = machine(io, X, y)
    fit!(m)
    p = predict(m, Xnew)
    @test p ≈ pred

end

end # module
true


text:
# http://www.voxforge.org/home/dev/acousticmodels/linux/create/htkjulius/tutorial/data-prep/step-2

###############################################################################
#
#    Copyright (C) 2015  VoxForge
#
#    This program is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
###############################################################################

function prompts2wlist(prompts, wlist)
  if ! isfile(prompts)
    error("can't find prompts file: $prompts")
  end

  wordhash = Dict{String, Int32}()
  prompts_fh=open(readlines, prompts)
  for lineln=prompts_fh
    line=chomp(lineln)
    line_array=split(line,r"\s+");
    shift!(line_array)
    for word=line_array
      wordhash[word]=1
    end
  end
  wordhash["SENT-END"]=1
  wordhash["SENT-START"]=1

  wordlist = keys(wordhash) # returns an iterator
  wlist_arr=Array(String,length(wordhash))
  i=1
  for word=wordlist
    wlist_arr[i] = word * "\n"
    i=i+1
  end
  sortedwlist_arr=sort(wlist_arr)

  wlist_fh=open(wlist,"w");
  write(wlist_fh,sortedwlist_arr);
  close(wlist_fh)
end

# if called from command line
if length(ARGS) > 0
  if ! isfile(ARGS[1])
    error("can't find prompts file: $ARGS[1]")
  end
  if length(ARGS) <= 2
    prompts2wlist(ARGS[1],ARGS[2] )
  else
    error("prompts2list: too many arguments for call from command line")
  end

end

max_stars_repo_name: Hyperparticle/lct-master

max_issues_repo_name: Hyperparticle/lct-master

max_forks_repo_name: Hyperparticle/lct-master


text:
"""
    @log var_name = val
    @log var_name expr

Log `var_name` to the `GLOBAL_LOG`. If `@log` is placed on an assignment for `var_name`, the variable
will also be created in the local scope. If `@log` is placed with a variable name before an expression,
the expression will run in the local scope without creating the variable in that scope.

## Example
```julia
function lorenz!(du, u, p, t)
    @log a = u[2]-u[1]
    @log b u[3]+a
    du[1] = p[1]*a
    du[2] = u[1]*(p[2]-u[3]) - u[2]
    du[3] = u[1]*u[2] - p[3]*u[3]
end
```
In this example, `a` is evaluated into the scope of the `lorenz!` function and is able to be used within
that scope. `b`, however is not evaluated into the `lorenz!` scope, so no `b` variable is created in
that scope. Both `a` and `b` will be logged.

"""
macro log(expr)
    return if expr isa Symbol
        quote
            if is_active(GLOBAL_LOG)
                local val = $(esc(expr))
                local var_name = $((expr,))[1]
                if !haskey(value_dict(GLOBAL_LOG), var_name)
                    setproperty!(GLOBAL_LOG, var_name, typeof(val)[])
                end
                push!(getindex(value_dict(GLOBAL_LOG), var_name), val)
            end
            nothing
        end
    elseif expr.head == :(=)
        if expr.args[1] isa Expr && expr.args[1].head == :tuple
            quote
                if is_active(GLOBAL_LOG)
                    local vals = $(esc(expr.args[2]))
                    local var_names = $(esc(expr.args[1].args))
                    for (var_name, val) in zip(var_names, vals)
                        if !haskey(value_dict(GLOBAL_LOG), var_name)
                            setproperty!(GLOBAL_LOG, var_name, typeof(val)[])
                        end
                        push!(getindex(value_dict(GLOBAL_LOG), var_name), val)
                    end
                end
                $(esc(expr))
            end
        else
            quote
                if is_active(GLOBAL_LOG)
                    local val = $(esc(expr.args[2]))
                    local var_name = $((expr.args[1],))[1]
                    if !haskey(value_dict(GLOBAL_LOG), var_name)
                        setproperty!(GLOBAL_LOG, var_name, typeof(val)[])
                    end
                    push!(getindex(value_dict(GLOBAL_LOG), var_name), val)
                end
                $(esc(expr))
            end
        end
    else
        :(error("Must be logging to some variable, use either `@log var_name val` or `@log var_name = val` forms"))
    end
end

macro log(var_name, expr)
    quote
        local val = $(esc(expr))
        if is_active(GLOBAL_LOG)
            local var_name = $((var_name,))[1]
            if !haskey(value_dict(GLOBAL_LOG), var_name)
                setproperty!(GLOBAL_LOG, var_name, typeof(val)[])
            end
            push!(getindex(value_dict(GLOBAL_LOG), var_name), val)
        end
        val
    end
end

max_stars_repo_name: jonniedie/SimulationLogs

max_issues_repo_name: jonniedie/SimulationLogs

max_forks_repo_name: jonniedie/SimulationLogs


text:
using Parameters
using ForwardDiff
using NLsolve
using ConstLab
using Devectorize

using Voigt
using Voigt.Unicode

const II = veye(6) ⊗ veye(6)
const Idev = eye(6,6) - 1/3 * II

@with_kw immutable CrystPlastMS
    n_ε_p::Vector{Float64}
    # One per slip system:
    n_κ::Vector{Float64}
    n_τ::Vector{Float64}
    n_μ::Vector{Float64}
end

function CrystPlastMS(nslip)
    n_ε_p = zeros(6)
    n_κ = zeros(nslip)
    n_τ = zeros(nslip)
    n_μ = zeros(nslip)
    CrystPlastMS(n_ε_p, n_κ, n_τ, n_μ)
end

@with_kw immutable CrystPlastMP
    E::Float64
    Ee::Matrix{Float64}
    ν::Float64
    σy::Float64
    n::Float64
    H::Float64
    q::Float64
    D::Float64
    tstar::Float64
    angles::Vector{Float64}
    sxm_sym::Vector{Vector{Float64}}
    nslip::Int
end

function CrystPlastMP(E, ν, σy, n, H, q, D, tstar, angles)
    sxm_sym = Vector{Vector{Float64}}()
    nslip = length(angles)
    for α = 1:nslip
        t = deg2rad(angles[α])
        s = [cos(t), sin(t), 0.0]
        m = [cos(t + pi/2), sin(t + pi/2), 0.0]
        sxm = voigtsym(s * m')
        push!(sxm_sym, sxm)
    end
    G = E / 2(1 + ν);
    K = E / 3(1 - 2ν);
    Ee = 2 * G * Idev+ K * II;
    Ee[4:6, 4:6] /= 2.0;
    Ee
    CrystPlastMP(E, Ee, ν, σy, n, H, q, D, tstar, angles, sxm_sym, nslip)
end

function stress(ε, dt, matpar::CrystPlastMP, matstat::CrystPlastMS)
    @unpack matstat: n_ε_p, n_κ, n_τ, n_μ
    @unpack matpar: Ee, ν, σy, sxm_sym, nslip

    σ_tr = Ee * (ε - n_ε_p)

    x0 = [σ_tr; n_κ; n_μ]

    R!(x, fx) = compute_residual!(fx, x, σ_tr, dt, matpar, matstat)
    res = nlsolve(R!, x0; iterations = 30, ftol = 1e-4, autodiff = true)
    if !NLsolve.converged(res)
        error("No convergence in material routine")
    end
    X = res.zero::Vector{Float64}
    ConstLab.@unpack_cp (σ, κ, μ) = X
    ε_p = copy(n_ε_p)
    τ = zeros(nslip)
    for α = 1:nslip
        τ[α] = σ : sxm_sym[α]
        ε_p += μ[α] * sxm_sym[α] * sign(n_τ[α])
    end

    return σ, CrystPlastMS(ε_p, κ, τ, μ)
end

function compute_residual!(fx, x, σ_tr, dt, matpar, matstat)
    @unpack matstat: n_ε_p, n_κ, n_τ, n_μ
    @unpack matpar: Ee, ν, σy, sxm_sym, nslip, H, q, tstar, D, n

    ConstLab.@unpack_cp (σ, κ, μ) = x

    R_Φ = zeros(eltype(σ), nslip)
    R_κ = zeros(eltype(σ), nslip)
    Φ = zeros(eltype(σ), nslip)
    τ = zeros(eltype(σ), nslip)

    for α=1:nslip
        R_κ[α] = κ[α] - n_κ[α]
        for β = 1:nslip
            R_κ -= H * μ[β] * (q + (1 - q) * Int(α == β))
        end
        τ[α] = σ : sxm_sym[α]
        Φ[α] = abs(τ[α]) - (κ[α] + σy)
        R_Φ[α] = μ[α] * tstar - dt * max(0, Φ[α]/D)^n
    end

    R_σ = σ - σ_tr
    ep = zeros(eltype(μ), 6)
    for α = 1:nslip
        ep += μ[α] .* sxm_sym[α] .* sign(τ[α])
    end

    R_σ += Ee * ep

    ConstLab.@pack_cp fx = (R_σ, R_κ, R_Φ)
    return
end


text:
function formatinfo(s::Spec, x::Irrational)
    if !isspecified(s.type) || s.type == 's'
        # string representation using String
        val = string(x)
        size, meta = formatinfo(Spec(s, align = default(s.align, ALIGN_RIGHT), type = 's'), val)
    else
        # numeric representation using BigFloat
        precision = default(s.precision, 6)
        val = BigFloat(x; precision = precision10to2(precision))
        size, meta = formatinfo(s, val)
    end
    return size, (val, meta)
end

function formatfield(data::Vector{UInt8}, p::Int, s::Spec, ::Irrational, meta)
    if !isspecified(s.type) || s.type == 's'
        s = Spec(s, align = default(s.align, ALIGN_RIGHT), type = 's')
    end
    return formatfield(data, p, s, meta[1], meta[2])
end

precision10to2(prec) = ceil(Int, (prec + 1) * log2(10))


text:
module D1_Interval_Test

using Compat
using Compat.Test
using IntervalArithmetic
using EAGOBranchBound

X = [Interval(1,5)]
model = BnBModel(X)
solver = BnBSolver()
set_to_default!(solver)
solver.Verbosity = "Normal"
solver.opt = [x -> x[1]^2]
solver.Lower_Prob = (X,k,p,opt,temp) -> [opt[1](X).lo,mid.(X),true,[]]
solver.Upper_Prob = (X,k,p,opt,temp) -> [opt[1](X).hi,mid.(X),true,[]]
a = solver.Lower_Prob(X,1,1,solver.opt,[])
b = solver.Upper_Prob(X,1,1,solver.opt,[])
solveBnB!(solver,model)

@test -1E-4 <= getobjval(model)-1 <= 1E-4

solver.iter_lim = true
solver.max_iter = 10
solver.Lower_Prob = (X,k,p,opt,temp) -> [opt[1](X).lo,mid.(X),false,[]]
model2 = BnBModel(X)
solveBnB!(solver,model2)
c1 = getobjval(model2)
@test c1 == Inf

solver.Lower_Prob = (X,k,p,opt,temp) -> [opt[1](X).lo,mid.(X),true,[]]
solver.Upper_Prob = (X,k,p,opt,temp) -> [opt[1](X).hi,mid.(X),false,[]]
model3 = BnBModel(X)
solveBnB!(solver,model3)
c2 = getobjval(model3)
@test c2 == Inf

solver.Lower_Prob = (X,k,p,opt,temp) -> [opt[1](X).lo,mid.(X),true,[]]
solver.Upper_Prob = (X,k,p,opt,temp) -> [opt[1](X).hi,mid.(X),true,[]]
function test_check(S,M,X1,X2)
    return (M.lbcnt == 2)
end
solver.Repeat_Chk = test_check
model5 = BnBModel(X)
solveBnB!(solver,model5)
c5 = getobjval(model5)
@test c5 == 1.5625

solver.Preprocess =  (feas_Pre,nsBox,UBDg,k_int,pos,opt,LBDn,UBDn,x,y) -> (false,nsBox)
solver.Lower_Prob = (X,k,p,opt,temp) -> [opt[1](X).lo,mid.(X),true,[]]
solver.Upper_Prob = (X,k,p,opt,temp) -> [((k==2) ? opt[1](X).hi : Inf),mid.(X),false,[]]
model4 = BnBModel(X)
solveBnB!(solver,model4)
c3 = getobjval(model4)
@test c3 == Inf
end


text:
struct SearchNode{T} end

const fname = ["has_child", "has_parent" , "nested"]

function Base.push!(t::AbstractDict, b::AbstractDict)
	for (x,y) in b
		push!(t, x => y)
	end
end

function make_loop(exprs::Vector )

	len = length(exprs)
	val = Array{Expr}(undef, length(exprs))

	for i in 1:len
		(methods, name, content)  = estrans(exprs[i])
		isa(name, Symbol) && (name = string(name))
		val[i] =  :( $methods => Dict( $name  => $content) )
	end

	esc(Expr( :call, :Dict, val... ))
end

function make_loop2(exprs::Vector )

	len = length(exprs)
	val = Array{Expr}(undef, length(exprs))

	for i in 1:len
		( methods , name, content)  = estrans(exprs[i])
		isa(name, Symbol) && (name = string(name))
		val[i] = methods == nothing  ? :($name => $(content)[$name]) : :($name  => $content)
	end

	esc(Expr( :call, :Dict, val... ))
end

function make_json(method::AbstractString, exprs::Vector )

	len = length(exprs)
	val = Array{Expr}(undef, length(exprs))

	for i in 1:len
		(methods, name, content)  = estrans(exprs[i])
		isa(name, Symbol) && (name = string(name))

		if methods == nothing
			val[i] = name in fname ? :($name => $(content)[$name]) : :( "bool" => Dict( $name => $(content)[$name]))
		else
			val[i] =  :( $methods => Dict( $name  => $content))
		end
	end
 	Expr( :call, :vcat, val... ) |> df -> esc(:(Dict($method => $df )))
end

function make_json(  exprs::Vector, type::AbstractString )

	len = length(exprs)
	val = Array{Expr}(undef, length(exprs))

	for i in 1:len
		( methods , name, content)  = estrans(exprs[i])
		isa(name, Symbol) && (name = string(name))

		if methods == nothing
			val[i] =  :( $name => $(content)[$name] )
		else
			val[i] =  :( $name  => $content)
		end
	end
	Expr( :call, :Dict, val... ) |> df -> esc(:(Dict($type => $df)))
end

function make_json( exprs::Vector )

	len = length(exprs)
	val = Expr[]
	query = Expr[]

	for i in exprs
		( methods , name, content)  = estrans(i)
		isa(name, Symbol) && (name = string(name))

		if methods == nothing
			push!(query, :($name => $(content)[$name]  ))
		else
			push!(val, :($name => $content ))
		end
	end

	vexp = Expr( :call, :Dict, val... )

	if ! isempty(query)
		qexp = Expr( :call, :Dict, query... ) |> df -> :(Dict("query" => Dict("bool" => $df)))
		esc(Expr( :call, :merge, vexp , qexp ))
	else
		esc(vexp)
	end

end

macro comm( expr...)
	return  make_loop(collect( expr ))
end

macro smi( expr...)
	return  make_loop2(collect( expr ))
end

macro fulltext( expr... )
    return  make_json(collect( expr ) ,"query")
end

macro filter( expr... )
    return  make_json( "filter", collect( expr ))
end

macro should( expr... )
    return  make_json( "should", collect( expr ))
end

macro must( expr... )
    return make_json( "must", collect( expr ))
end

macro must_not( expr... )
    return  make_json( "must_not", collect( expr ))
end

macro query( expr... )
    return  make_json( collect( expr ))
end

macro nested( expr... )
    return  make_json(collect( expr ) ,"nested")
end

macro has_child( expr... )
    return  make_json(collect( expr ) ,"has_child")
end

macro has_parent( expr... )
    return  make_json(collect( expr ) ,"has_parent")
end

function estrans( expr::Expr )
    estrans( SearchNode{expr.head}, expr)
end

function estrans(expr::Symbol)
    exp =  quote
		try
			$expr.keys[$expr.slots .!= 0][1]
		catch
			string($(Expr(:quote, expr)))
		end
	   end
    ("Symbol"  , exp , expr )
end

function estrans(::Type{SearchNode{:(*)}}, expr::Expr)
	("regexp" , expr.args[2] , expr.args[3] )
end

function estrans(::Type{SearchNode{:(%)}}, expr::Expr)
	("wildcard" , expr.args[2] , expr.args[3] )
end

function estrans(::Type{SearchNode{:(|)}}, expr::Expr)
	("prefix" , expr.args[2] , expr.args[3] )
end

function estrans( ::Type{SearchNode{:call}}, expr::Expr)
     estrans(SearchNode{expr.args[1]}, expr)
end

function estrans(::Type{SearchNode{:(=)}}, expr::Expr)
	("term" , expr.args[1] ,  expr.args[2] )
end

function estrans(::Type{SearchNode{:(in)}}, expr::Expr)
	("terms" , expr.args[2] , expr.args[3] )
end

function estrans(::Type{SearchNode{:(has)}}, expr::Expr)
	("exists" , "field" , expr.args[2] )
end

function estrans(::Type{SearchNode{:(>=)}}, expr::Expr)
	("range" , expr.args[2] , :(Dict( "gte" => $(expr.args[3])) ))
end

function estrans(::Type{SearchNode{:(<=)}}, expr::Expr)
	("range" , expr.args[2] , :(Dict( "lte" => $(expr.args[3])) ))
end

function estrans(::Type{SearchNode{:(>)}}, expr::Expr)
	("range" , expr.args[2] , :(Dict( "gt" => $(expr.args[3])) ))
end

function estrans(::Type{SearchNode{:(<)}}, expr::Expr)
	("range" , expr.args[2] , :(Dict( "lt" => $(expr.args[3])) ))
end

function estrans(::Type{SearchNode{:comparison}}, expr::Expr)
	sym = ( syms(esyms{expr.args[2]}) ,syms( ! esyms{expr.args[4]}) )
	("range" , expr.args[3] , :(Dict( $(sym[1]) => $(expr.args[1]) ,$(sym[2]) => $(expr.args[5]) )))
end

abstract type esyms{T} end

syms(::Type{esyms{:<}})       = "gt"
syms(::Type{esyms{:(<=)}})    = "gte"
syms(::Type{esyms{:>}})       = "lt"
syms(::Type{esyms{:(>=)}})    = "lte"

Base.:!(::Type{esyms{:<}})    = esyms{:>}
Base.:!(::Type{esyms{:(<=)}}) = esyms{:(>=)}
Base.:!(::Type{esyms{:>}})    = esyms{:<}
Base.:!(::Type{esyms{:(>=)}}) = esyms{:(<=)}

function estrans(::Type{SearchNode{:macrocall}}, expr::Expr)
	(nothing  , replace(string(expr.args[1]), "@" => "")  , expr)
end

max_stars_repo_name: frankill/ES

max_issues_repo_name: frankill/ES

max_forks_repo_name: frankill/ES


text:
################################################################################
#
#    NfOrd/TorsionUnits.jl : Torsion units in generic number field orders
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016 Tommy Hofmann
#
################################################################################

export torsion_unit_group

################################################################################
#
#  Torsion unit test
#
################################################################################

@doc Markdown.doc"""
    istorsion_unit(x::NfOrdElem, checkisunit::Bool = false) -> Bool

Returns whether $x$ is a torsion unit, that is, whether there exists $n$ such
that $x^n = 1$.

If `checkisunit` is `true`, it is first checked whether $x$ is a unit of the
maximal order of the number field $x$ is lying in.
"""
function istorsion_unit(x::NfOrdElem, checkisunit::Bool = false)
  return istorsion_unit(x.elem_in_nf, checkisunit)
end

################################################################################
#
#  Order of a single torsion unit
#
################################################################################

@doc Markdown.doc"""
    torsion_unit_order(x::NfOrdElem, n::Int)

Given a torsion unit $x$ together with a multiple $n$ of its order, compute
the order of $x$, that is, the smallest $k \in \mathbb Z_{\geq 1}$ such
that $x^k = 1$.

It is not checked whether $x$ is a torsion unit.
"""
function torsion_unit_order(x::NfOrdElem, n::Int)
  return torsion_unit_order(x.elem_in_nf, n)
end

################################################################################
#
#  Interface: number fields
#
################################################################################

function istorsion_unit_group_known(K::AnticNumberField)
  return _get_nf_torsion_units(K, false) != nothing
end

function torsion_unit_group(K::AnticNumberField)
  ord, gen = _torsion_units_gen(K)
  mp = AbToNfMultGrp(K, ord, gen)
  return domain(mp), mp
end

function torsion_units(K::AnticNumberField)
  ord, gen = _torsion_units_gen(K)
  return powers(gen, ord-1)
end

function torsion_units_generator(K::AnticNumberField)
  ord, gen = _torsion_units_gen(K)
  return gen
end

function torsion_units_order(K::AnticNumberField)
  ord, gen = _torsion_units_gen(K)
  return ord
end

function torsion_units_gen_order(K::AnticNumberField)
  ord, gen = _torsion_units_gen(K)
  return gen, ord
end

################################################################################
#
#  Interface: orders
#
################################################################################

@doc Markdown.doc"""
    torsion_units(O::NfOrd) -> Array{NfOrdElem, 1}

Given an order $O$, compute the torsion units of $O$.
"""
function torsion_units(O::NfOrd)
  g, ord = torsion_units_gen_order(O)
  return powers(g, ord-1)
end

@doc Markdown.doc"""
    torsion_units_generator(O::NfOrd) -> NfOrdElem

Given an order $O$, compute a generator of the torsion units of $O$.
"""
function torsion_units_generator(O::NfOrd)
  return torsion_units_gen_order(O::NfOrd)[1]
end

@doc Markdown.doc"""
    torsion_units_gen_order(O::NfOrd) -> NfOrdElem

Given an order $O$, compute a generator of the torsion units of $O$ as well as its order.
"""
function torsion_units_gen_order(O::NfOrd)
  ord, g = _torsion_units_gen(nf(O))
  if ismaximal_known_and_maximal(O)
    return O(g), ord
  end
  #We need to check which torsion units are in the order.
  lf = factor(ord)
  ord_O = 1
  for (p, v) in lf
    if p == 2 && isone(v)
      ord_O *= 2
      continue
    end
    p_int = Int(p)
    ord_p = p_int^v
    u = g^divexact(ord, ord_p)
    while !(u in O)
      u = u^p_int
      ord_p = divexact(ord_p, p_int)
    end
    ord_O *= ord_p
  end
  #ord_O is the order of the torsion units in O
  @assert iszero(mod(ord, ord_O))
  genO = O(g^divexact(ord, ord_O))
  return genO, ord_O
end

@doc Markdown.doc"""
    torsion_unit_group(O::NfOrd) -> GrpAb, Map

Given an order $O$, returns the torsion units as an abelian group $G$
together with a map $G \to \mathcal O^\times$.
"""
function torsion_unit_group(O::NfOrd)
  g, ord = torsion_units_gen_order(O)
  f = AbToNfOrdMultGrp(O, ord, O(g))
  return domain(f), f
end

################################################################################
#
#  Torsion units via lattice enumeration
#
################################################################################

function _torsion_units_lattice_enum(O::NfOrd)
  n = degree(O)
  K = nf(O)
  r1, r2 = signature(K)
  B = basis(O)

  if r1 > 0
    return [ O(1), -O(1) ], -O(1)
  end

  function _t2_pairing(x, y, p)
    local i
    v = minkowski_map(x, p)
    w = minkowski_map(y, p)

    t = zero(parent(v[1]))

    for i in 1:r1
      t = t + v[i]*w[i]
    end

    for i in (r1 + 1):(r1 + 2*r2)
      t = t + v[i]*w[i]
    end

    return t
  end

  p = 64

  gram_found = false

  could_enumerate = false

  A = ArbField(p, false)
  M = ArbMatSpace(A, n, n,false)()

  while true
    A = ArbField(p, false)
    M = ArbMatSpace(A, n, n, false)()

    gram_found = true

    @vprint :UnitGroup 1 "Computing Gram matrix with precision $p ... \n"

    for i in 1:n, j in 1:n
      M[i,j] = _t2_pairing(B[i], B[j], p)
      if !isfinite(M[i, j])
        p = 2*p
        gram_found = false
        break
      end
    end

    if gram_found
      break
    end
  end

  @vprint :UnitGroup 1 "Gram matrix is $M\n"

  @vprint :UnitGroup 1 "Enumerating elements with T_2 bounded by $n ... \n"
  l = enumerate_using_gram(M, A(n))

  R = Array{NfOrdElem, 1}()

  for i in l
    if O(i) == zero(O)
      continue
    end
    if istorsion_unit(O(i))
      push!(R, O(i))
    end
  end

  i = 0

  for outer i in 1:length(R)
    if torsion_unit_order(R[i], length(R)) == length(R)
      break
    end
  end

  return R, R[i]
end

################################################################################
#
#  Torsion units via root lifting
#
################################################################################

# For each n _euler_phi_inverse_maximum_ is the maximal i,
# such that varphi(i) divides n
#
# Magma code:
# [Maximum(&cat[ EulerPhiInverse(d) : d in Divisors(n)  | #EulerPhiInverse(d) ne 0 ]) : n in [1..250]];
const _euler_phi_inverse_maximum =
[ 2, 6, 2, 12, 2, 18, 2, 30, 2, 22, 2, 42, 2, 6, 2, 60, 2, 54, 2, 66, 2, 46, 2,
90, 2, 6, 2, 58, 2, 62, 2, 120, 2, 6, 2, 126, 2, 6, 2, 150, 2, 98, 2, 138, 2,
94, 2, 210, 2, 22, 2, 106, 2, 162, 2, 174, 2, 118, 2, 198, 2, 6, 2, 240, 2,
134, 2, 12, 2, 142, 2, 270, 2, 6, 2, 12, 2, 158, 2, 330, 2, 166, 2, 294, 2, 6,
2, 276, 2, 62, 2, 282, 2, 6, 2, 420, 2, 6, 2, 250, 2, 206, 2, 318, 2, 214, 2,
378, 2, 242, 2, 348, 2, 18, 2, 354, 2, 6, 2, 462, 2, 6, 2, 12, 2, 254, 2, 510,
2, 262, 2, 414, 2, 6, 2, 274, 2, 278, 2, 426, 2, 6, 2, 630, 2, 6, 2, 298, 2,
302, 2, 30, 2, 46, 2, 474, 2, 6, 2, 660, 2, 486, 2, 498, 2, 334, 2, 588, 2, 22,
2, 346, 2, 118, 2, 690, 2, 358, 2, 594, 2, 6, 2, 564, 2, 18, 2, 12, 2, 382, 2,
840, 2, 6, 2, 394, 2, 398, 2, 750, 2, 6, 2, 618, 2, 6, 2, 636, 2, 422, 2, 642,
2, 6, 2, 810, 2, 6, 2, 726, 2, 446, 2, 870, 2, 454, 2, 458, 2, 94, 2, 708, 2,
158, 2, 12, 2, 478, 2, 1050, 2, 46, 2, 12, 2, 166, 2, 30, 2, 502 ]

# One should/could also try to be closer to Algorithm 1
# in Molin, "On the calculation of roots of unity in a number field"
function _torsion_group_order_divisor(K::AnticNumberField)

  if degree(K) <= 250
    upper_bound = _euler_phi_inverse_maximum[degree(K)]
  else
    error("Not implemented yet")
  end

  p = upper_bound + 1
  m = fmpz(0)
  stable = 0

  first = true

  if ismaximal_order_known(K)
    disc = abs(discriminant(maximal_order(K)))
  elseif isdefining_polynomial_nice(K)
    disc = discriminant(EquationOrder(K))
  else
    disc_1 = discriminant(K.pol)
    disc = numerator(disc_1)*denominator(disc_1)
  end
  threshold = 5*degree(K)

  while true
    p = next_prime(p)
    Rp = Nemo.GF(p, cached=false)
    Rpt, t = PolynomialRing(Rp, "t", cached=false)
    gp = Rpt(K.pol)

    if degree(gp) != degree(K) || !issquarefree(gp)
      continue
    end

    lp = collect(keys(factor_distinct_deg(gp)))

    minf = minimum(lp)

    if first
      m_new = fmpz(p)^minf - 1
      m_new, _ = ppio(m_new, disc)
      if isodd(m_new)
        m_new = 2 * m_new
      end
      first = false
    else
      m_new =  gcd(m, powermod(fmpz(p), minf, m) - 1)
    end

    if m_new == 2
      return fmpz(2)
    end

    if m_new == m
      stable += 1
    else
      stable = 0
    end
    if !divisible(fmpz(degree(K)), euler_phi(m_new))
      stable = 0
    end

    if stable == threshold && m_new <= upper_bound
      return m_new
    end

    m = m_new
  end
end

function _torsion_units_gen(K::AnticNumberField)
  if istorsion_unit_group_known(K)
    c = _get_nf_torsion_units(K)::Tuple{Int, nf_elem}
    return c[1], c[2]
  end
  r1, r2 = signature(K)
  if r1 > 0
    return 2, K(-1)
  end

  m = _torsion_group_order_divisor(K)
  Ky = PolynomialRing(K, "y", cached = false)[1]
  fac = factor(m).fac
  gen = K(1)
  ord = 1
  Zx, x = PolynomialRing(FlintZZ, "x")
  for (p, v) in fac
    if p == 2 && v == 1
      mul!(gen, gen, K(-1))
      ord *= 2
      continue
    end
    for i = v:-1:1
      f = cyclotomic(Int(p)^i, x)
      fK = map_coefficients(K, f, parent = Ky)
      r = _roots_hensel(fK, max_roots = 1, isnormal = true, root_bound = fmpz[one(fmpz) for i in 1:(r1 + r2)])
      if length(r) > 0
        mul!(gen, gen, r[1])
        ord *= Int(p)^(i)
        break
      end
    end
  end
  _set_nf_torsion_units(K, (ord, gen))
  return ord, gen
end