Julia 101 – 0.1. Least-Square fitting

year = Float64[0]

for i in 1:6
  push!(year, i)

year = year .+ 2013
7-element Array{Float64,1}:
year_x = collect(Float64, 2013:1:2020);
h = [2, 2, 2, 3, 4, 6, 8, 12];
using Plots
scatter(year_x, h, 
  xlab = "year", ylab = "h-index", 
  xlim = [2012, 2023], ylim = [0,15], framestyle = :box)

scatter!(year_x[3:end], h[3:end])
using LsqFit

@. model(x, p) = p[1] + x*p[2] + x^2.0*p[3]  
fit = curve_fit(model, year_x, h, [2.0, 0.5, -3.0])
LsqFit.LsqFitResult{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Int64,1}}([1.1832864618761656e6, -1174.946148708822, 0.2916665972735931], [0.20833284547552466, -0.19642864260822535, -0.017856936203315854, -0.2559520348440856, 0.08928606053814292, 0.017857350641861558, 0.5297618352342397, -0.3750004852190614], [1.0000000000082452 2012.9999999932165 4.05216899997484e6; 1.0000000000082452 2014.0000000073403 4.05619599998367e6; … ; 1.0000000000082452 2018.9999999797847 4.076360999975335e6; 1.0000000000082452 2019.9999999939084 4.08039999999393e6], true, Int64[])
T = typeof(fit)
for (name, typ) in zip(fieldnames(T), T.types)
    println("type of the fieldname name istyp")
type of the fieldname param is Array{Float64,1}
type of the fieldname resid is Array{Float64,1}
type of the fieldname jacobian is Array{Float64,2}
type of the fieldname converged is Bool
type of the fieldname wt is Array{Int64,1}
search: fit LsqFit filter filter! curve_fit first firstindex isfinite popfirst!

No documentation found.

fit is of type LsqFit.LsqFitResult{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Int64,1}}.


struct LsqFit.LsqFitResult{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Int64,1}} <: Any


param     :: Array{Float64,1}
resid     :: Array{Float64,1}
jacobian  :: Array{Float64,2}
converged :: Bool
wt        :: Array{Int64,1}
fitted_p = fit.param
3-element Array{Float64,1}:
new_x = collect(Float64, 2013:1:2025)
new_y = fitted_p[1] .+ new_x.*fitted_p[2] .+ (new_x.^2).*fitted_p[3] 
plot!(new_x, new_y, ylim = [0,30], seriestype =:line)
8-element Array{Float64,1}:
8×3 Array{Float64,2}:
 1.0  2013.0  4.05217e6
 1.0  2014.0  4.0562e6
 1.0  2015.0  4.06022e6
 1.0  2016.0  4.06426e6
 1.0  2017.0  4.06829e6
 1.0  2018.0  4.07232e6
 1.0  2019.0  4.07636e6
 1.0  2020.0  4.0804e6
0-element Array{Int64,1}
Astronomy · computer · Linux · Note


有时候你想在一个目录的所有子目录下批量建立文件夹,并且批量的运行这些子目录内新建立好的文件夹内的 IDL 程序,以下是一种可能的解决途径,如果你有更好的方法,欢迎告诉我

Here are the things I want to do:

  1.  I have a data folder containing many sources, like Arp220, NGC981, etc.. Each of them possesses a folder under the whole data folder (./source/Arp220, ./source/NGC891, …). I want to create a subfolder named “plot” under each sources folder, e.g., ./source/Arp220/plot.
  2. And then I want to copy the IDL routine to each “plot” folder I just created and run all the PRO routines at one time.

Here is a bash script being able to do such things:

for i in $(ls -d */); do mkdir -p $i/'plot'; done
for i in $(ls -d */); do cp ./ $i/'plot'/; done
for i in $(ls -d */); do cd $PWD/$i/'plot';idl -e ".r ./"; cd ../.. ; done

The first line means creating a “plot” folder under all the subfolders under you present folder. And the next line is to copy the IDL routine ( to each “plot” folder you created. The last line runs all the IDL routines.


find ./* -type f|grep X|while read file; do rm "$file"; done

比如你要删除“./you suck/arXiv.pdf“这种文件。上面那句话中 -type f|grep X|while read file 是为了保证有空格的目录也可以被写入到一个变量里面


find ./* -type f -name "*X*" -exec rm {} \;
Astronomy · computer · Internet · Linux · Note

Obfuscated ssh Client under Linux

SSH is a powerful tool not only for astronomers. We often use it to download observation data from the observatory’s server or run some complicated programs on remote (overseas) servers. But for the reason as we all know, the SSH tunnels are always being obstructed. So we need some technique to get over such obstruction. Therefore, the obfuscated SSH turns out to be a sufficient method.

I will describe how to use obfuscated SSH_TUNNEL under Linux as below:
(Original Code: Here)


cd ./brl-obfuscated*
./configure  --prefix=/usr/local/newssh --sysconfdir=/etc/newssh
make install

Then you need to configure the ssh server following the guide here.

Then use this to run the ssh client:

/usr/local/newssh/bin/ssh -N -v -Z ObfuscateKeyWord -p ***** username@hostname -D

There are 4 parameters you need to specify: “ObfuscateKeyWord” is the obfuscated keyword you set in the configuration file (sshd_config) of the server, “*****” is the ssh server port number, “username” and “hostname” are as the name described.

Astronomy · computer · Note

Least-Squares Fitting of data containing errors

(This is a note personal only! )

There is a powerful tool dealing with data containing both “X” and “Y” errors to make a lease-square fitting in IDL, that is mpfit. This routine allows you to use a self-defined function to do the fitting. Also, you can set up some conditions such as the parameters’ limit or fixed parameters etc. This is fabulous! So beginning with the actual astronomical data have errors both in “X” and “Y”, you need to use a function called “linfitex” (you can see more info from the mpfit website) to both define the function you want to fit and calculate the least-square needed in the fitting. Here is an example of the function:

FUCNTION myfunct, p, x=x, y=y \$
     , sigma_x=sigma_x, sigma_y=sigma_y, _EXTRA=extra
  a = p[0]
  b = p[1]
  f = a*alog10(X)+b
;#########calculate the derivative of the function#########
;##put into resid as (y-f)/sqrt(sigma_y^2+(y'*sigma_x)^2)## 
  resid = (y - f)/sqrt(sigma_y^2 + (a/(x*alog(10))*sigma_x)^2) 
RETURN, resid

In the code above I just defined a function as

    \[y=a \log_{10} (b)+c\]

Then you can read your data. and when you want to fit the data, just using mpfit as:

p      = [1,1,1]
p[1]   = 1         
p[2]   = 2
p[3]   = 3 ;initial guess of the parameters
param  = MPFIT('MYFUNCT', params,\$
      FUNCTARGS={X:x,Y:y,SIGMA_X:x_er,SIGMA_Y:y_er}, PERROR=p_er)

The FUNCTARGS is the place where your data go into the function. And then you can find the results of the parameters from the fitting by print “param” and also the 1σ error from printing p_er. If you only have the error on “Y”, you can use a function like this:

FUNCTION myfunct, p, x=x, y=y, er=er
  a = p[0]     
  b = p[1]
  c = p[2]
  f = a*alog10(b)+c
RETURN, (y-f)/er

And then using the same way described above to call “mpfit” and get the results and errors.

Here is also a PYTHON version of this procedure.