Julia 101 – 0.1. Least-Square fitting

July 26, 2020 by starwave

year = Float64[0]

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

year = year .+ 2013

7-element Array{Float64,1}:
 2013.0
 2014.0
 2015.0
 2016.0
 2017.0
 2018.0
 2019.0

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 typ")
end

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}

?fit

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}}.

Summary

struct LsqFit.LsqFitResult{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Int64,1}} &lt;: Any

Fields

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}:
     1.1832864618761656e6
 -1174.946148708822
     0.2916665972735931

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)

fit.resid

8-element Array{Float64,1}:
  0.20833284547552466
 -0.19642864260822535
 -0.017856936203315854
 -0.2559520348440856
  0.08928606053814292
  0.017857350641861558
  0.5297618352342397
 -0.3750004852190614

fit.jacobian

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

fit.wt

0-element Array{Int64,1}

fit.converged

true

利用BASH脚本运行IDL程序、创建子目录或删除文件

June 5, 2013June 27, 2017 by starwave

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

Here are the things I want to do:

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.
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:

#!/bin/bash
for i in $(ls -d */); do mkdir -p $i/'plot'; done
for i in $(ls -d */); do cp ./r_name.pro $i/'plot'/; done
for i in $(ls -d */); do cd $PWD/$i/'plot';idl -e ".r ./r_name.pro"; 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 (r_name.pro) to each “plot” folder you created. The last line runs all the IDL routines.

顺便再提一个操作，即“批量删除某个目录下所有子目录中包含”X”关键词的文件”。这个可以用一句话可以解决：

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

March 10, 2013June 27, 2017 by starwave

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)

INSTALLATION:

wget http://aenes.googlecode.com/files/brl-obfuscated-openssh.zip
unzip brl-obfuscated-openssh.zip
cd ./brl-obfuscated*
./configure  --prefix=/usr/local/newssh --sysconfdir=/etc/newssh
make
make install

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

USAGE:
Then use this to run the ssh client:

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

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

March 19, 2012June 27, 2017 by starwave

(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
END

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.

Rhapsody of the Universe

Category: computer